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1.0  INTRODUCTION 


Computer  aided  design  (CAD)  has  become  a  standard  tool  used  both  by  the  Navy 
and  contractors  in  the  design  of  ships.  Detailed  ship  drawings  are  generated  with  CAD 
software  such  as  AutoCAD.  The  advantage  is  that  modification  and  redesign  are  rel¬ 
atively  easy  because  of  the  associativity  of  entities  (i.e.,  points,  lines,  circles,  surfaces, 
etc.).  The  CAD  database  contains  the  relationships  between  all  entities.  Therefore, 
the  secondary  effects  of  modifying  an  edge,  for  example,  are  automatically  accounted 
for  when  other  structures  that  contain  the  modified  edge  are  generated. 

The  ability  for  rapidly  modifying  a  structure  is  ideal  for  concurrent  engineering, 
that  is,  the  simultaneous  design  of  a  system  by  all  engineering  disciplines.  Electro¬ 
magnetic  (EM),  fluid,  and  structural  designers  can  all  use  the  same  CAD  database 
for  their  analyses.  Any  changes  made  to  the  system  are  immediately  available  to  all 
engineers. 

For  the  electromagnetic  design  of  systems,  several  computer  codes  are  available 
to  government  agencies  and  contractors  that  use  triangular  facet  models  of  bodies 
to  define  scattering  and  radiation  structures.  They  include  PATCH,  CARLOS-3D, 
and  FERM  [1-3],  all  of  which  are  written  in  the  FORTRAN  language.  FERM  has  its 
own  preprocessing  program  to  generate  geometry  files.  Furthermore,  it  uses  binary 
data  files  which  are  not  easily  transportable  between  computer  codes  and  systems. 
CARLOS-3D  can  use  geometry  files  generated  by  the  CAD  program  ACAD  without 
modification.  However,  CARL0S-3D  can  only  solve  scattering  problems,  not  radia¬ 
tion  (antenna)  problems.  PATCH  is  capable  of  solving  both  scattering  and  radiation 
problems.  It  uses  an  ASCII  input  file  that  defines  the  body  on  the  basis  of  edge 
connections  (as  opposed  to  triangle  connections).  The  edge  definition  approach  can 
be  ambiguous  in  some  isolated  (but  predictable)  special  cases. 

For  the  primary  application  considered  here,  high  frequency  (HF)  ship  antenna 
analysis,  PATCH  has  been  found  to  be  most  useful.  However,  PATCH  has  only  a 
basic  geometry  preprocessor  capable  of  simple  shapes  such  as  plates,  cylinders,  cones, 
and  spheres.  The  data  for  more  complex  shapes  must  be  input  by  hand.  It  involves 
defining  the  nodes  of  every  triangle  on  the  body  and  specifying  an  edge  connection 
list.  Ship  models  that  are  meshed  for  use  up  to  30  MHz  have  about  3000  nodes 
(triangle  vertices)  and  5000  edges.  In  this  case,  data  entry  itself  is  a  major  effort. 
Furthermore,  if  a  ship  modification  is  to  be  investigated,  portions  of  the  structure 
may  have  to  be  remeshed  and  re- input,  which  again  represents  a  major  effort.  It 
is  apparent  that  for  computational  EM  codes  to  be  practical  engineering  tools,  an 
efficient  surface  meshing  operation  must  be  available. 
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The  CAD  application  ACAD  (Advanced  Computer  Aided  Design)  is  capable  of 
performing  the  desired  automeshing.  Databases  from  other  CAD  programs  can  be 
imported  into  ACAD  and  then  meshed  and  output  in  a  special  “facet”  format.  The 
facet  file  is  ASCII  and  contains  the  node  and  facet  information  required  by  PATCH, 
although  it  is  not  is  the  proper  format.  An  ACAD-to-PATCH  translator  was  written 
to  reformat  the  facet  file  into  one  that  is  recognized  by  PATCH.  Therefore,  it  is  now 
possible  to  take  a  NAVSEA  CAD  file,  mesh  it  automatically  in  ACAD,  and  then  use 
the  output  to  run  PATCH. 

This  report  summarizes  the  development  of  the  automeshing  process  and  the 
step-by-step  procedure  to  go  from  a  CAD  drawing  file  to  PATCH  input  file.  First  an 
overview  is  given  of  the  two  CAD  applications  used  here:  ACAD  and  AutoCAD.  The 
automeshing  procedure  is  described,  and  special  file  format  translation  software  is 
also  presented.  Finally,  some  shipboard  applications  are  discussed  and  a  few  helpful 
hints  and  guidelines  are  given. 

2.0  CAD  PROGRAMS 

2.1  ACAD  GENERAL  DESCRIPTION1 

ACAD  (Advanced  Computer  Aided  Design)  [4-5]  provides  users  with  the  abil¬ 
ity  to  create  and  modify  geometry  in  two-  or  three-dimensions.  Users  can  choose 
to  model  geometry  with  wireframes,  surfaces,  or  solids.  ACAD  is  the  primary  tool 
used  by  Lockheed  Fort  Worth  Company’s  Advanced  Programs  for  configuration  and 
subsystem  design  of  new  and  existing  aircraft  programs.  ACAD’s  primary  role  is  the 
generation  of  geometry  and  some  limited  analysis.  Much  of  the  analysis  performed 
within  ACAD  is  geometrical  analysis.  For  other  types  of  analysis,  ACAD  generates 
interface  files  for  transferring  to  groups  who  specialize  in  a  particular  analysis  field 
such  as  Radar  Cross  Section  (RCS),  Aero,  or  Computational  Fluid  Dynamics  (CFD). 

Inputting  data  to  ACAD  is  accomplished  through  one  of  many  input  modes  avail¬ 
able  to  the  designer.  Example  options  include  digitizing  locations,  entering  explicit 
coordinate  values,  snap  to  grid,  and  intersections.  Each  entity  (splines,  lines,  points, 
surfaces,  etc.)  can  have  individual  color,  width,  and  style  attributes.  Logical  group¬ 
ings  of  entities  can  be  separated  and  managed  with  layers,  groups,  and  blanking. 
ACAD  models  can  be  viewed  orthographically  or  in  perspective.  Users  can  specify 
view  orientation  and  choose  to  display  geometry  in  multiple  window  configurations. 
Window  operations  such  as  panning,  zooming,  and  auto  extents  are  accomplished  at 
any  time  providing  instream  capability.  The  ACAD  user  can  also  control  the  display 

1Most  of  the  material  in  this  section  is  taken  directly  from  reference  4 
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of  surfaces  or  solids  with  options  as  wireframe,  hidden  line  removal,  flat,  or  Gouraud 
shading. 

At  the  heart  of  the  ACAD  system  is  the  associative  database.  In  an  associa¬ 
tive  database,  geometry  is  linked  together  in  a  relational  structure  that  remembers 
parent/ child  dependencies.  This  type  of  database  enables  rapid  modifications  of  ge¬ 
ometry,  since  modifying  one  geometric  element  automatically  adjusts  its  dependencies 
based  on  a  set  of  predefined  rules.  For  instance,  changing  a  control  spline  of  a  fuse¬ 
lage  will  automatically  regenerate  any  surface(s)  built  with  the  spline.  In  turn,  any 
geometry  that  is  associated  to  the  fuselage  surface  (i.e.,  plane/curve  and  surface  in¬ 
tersections,  fillets)  will  automatically  regenerate. 

Within  ACAD  exists  a  read/ write  Initial  Graphics  Exchange  Specification  (IGES) 
translator.  The  IGES  translator  allows  ACAD  the  ability  to  exchange  drawing  data 
with  other  CAD  systems  that  support  IGES.  Example  CAD  systems  include  CAD- 
CAM,  CATIA,  COMPUTERVISION,  and  AutoCAD. 

The  ACAD  system  also  supports  a  binary  data  file  converter.  This  built  in  con¬ 
verter  enables  a  binary  file  created  on  the  SUN  to  be  read  in  directly  on  a  Silicon 
Graphics  or  Apollo  workstation  without  having  to  convert  to  a  neutral  ASCII  file. 
This  utility  is  extremely  beneficial  to  projects  supporting  a  mixture  of  workstations 
and  using  a  transparent  networking  system.  Hardcopy  output  is  available  on  a  variety 
of  devices  accessed  through  ACAD.  Additional  features  of  the  system  are: 

1.  Three-dimensional  lines  drawings 

2.  Analysis  geometry  models 

3.  Three-  and  Five- View  drawings 

Table  1  contains  more  information  on  the  types  of  commands  available  on  ACAD 
Version  9.0. 

2.2  AUTOCAD  GENERAL  DESCRIPTION  2 

The  capabilities  of  AutoCAD3  [6]  are  similar  to  those  of  ACAD.  AutoCAD  is  one 
of  the  most  widely  distributed  CAD  programs,  and  runs  on  all  platforms  (UNIX, 
DOS  and  Windows,  and  Macintosh).  As  in  the  case  of  ACAD,  the  effect  of  every 
change  made  to  a  drawing  appears  immediately  on  screen. 


2  Most  of  the  material  in  this  section  is  taken  directly  from  reference  6 

3 Frequently  the  AutoCAD  suite  of  programs  is  collectively  referred  to  as  ACAD.  This  terminology 
is  not  used  here  to  avoid  confusion  with  the  ACAD  discussed  in  Section  2.1. 
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Table  1:  Summary  of  ACAD  Commands 


Transformations 

Scale 

Translate 

Rotate 

Mirror 

Copy 

Display  Options 

1026  User  Defined  Layers 
Blank  On,  Off,  and  On  Only 
Color,  Style,  and  Width  Line  Fonts 
Hidden  Line,  Flat,  and  Gouraud  Shading 
Auxiliary  Viewing 

Orthographic  or  Perspective  Viewing 
Multiple  Windows  (up  to  6) 

Dynamic  Viewing 

Zooming,  Panning,  and  Auto  Extents 

Input  Options 

Digitize 

Reference  Existing  Data  Points 

Key  in  Explicit  x,  y,  z 

Intersection 

Point  On 

Snap  to  Grid 

Hierarchical  Input  Mode 


Intersections 
Curve-Curve 
Curve-Plane 
Curve- Surface 
Plane-Surface 
Surface-Surface 

Curve  Projections  onto  Surfaces 

Drafting  Utilities 

Break,  Trim,  and  Join  Curves 

Corner 

Grouping 

Construction  Planes 
Local  Coordinate  Systems 
Offsets 

Text  and  Dimensions 
Groups,  Dittos,  and  Details 
Crosshatching 

Three  Dimension  Design 

Point,  Line,  and  Spline  Primitives 
Conic,  Circles,  and  Ellipses  Primitives 
Six  Forms  of  Surfaces 
Curve  and  Surface  Editing 
Trimmed  Surfaces  (Faces) 

Mass  Properties  (volumes,  CGS,  areas) 
Offset  Surfaces 

Wireframe,  Surface,  k  Solid  Modeling 


< 
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AutoCAD  functions  lets  the  user  modify  the  drawing  in  a  variety  of  ways.  En¬ 
tities  can  be  erased  or  moved,  or  copied  to  form  repeated  patterns.  The  user  can 
change  the  view  of  the  drawing  displayed  on  screen,  or  display  information  about  the 
drawing.  AutoCAD  also  provides  drawing  aids  that  allow  the  positioning  of  entities 
accurately.  The  simple  command  format  of  AutoCAD  allows  the  user  to  accomplish 
most  of  the  functions  in  Table  1.  One  important  exception  is  the  inability  to  generate 
a  shell  mesh.  Therefore,  automeshing  must  be  performed  in  ACAD. 

2.3  IGES  FILES 

2.3.1  INTRODUCTION 

The  AutoCAD  and  ACAD  support  translation  of  drawings  to  and  from  the  Initial 
Graphics  Exchange  Specification  (IGES).  IGES  is  a  public-domain  data  specification 
intended  as  an  international  standard  for  the  exchange  of  information  between  CAD 
systems.4  An  IGES  translator  is  written  specifically  for  a  given  CAD  program  but, 
in  principle,  enables  drawings  to  be  transferred  to  and  from  other  CAD  systems  that 
also  support  IGES  (in  IGES,  these  are  known  as  sending  and  receiving  systems). 

Files  created  using  the  IGES  format  will  consist  of  at  most  six  sections  of  in¬ 
formation,  five  of  which  are  mandatory.  These  sections  must  exist  in  the  following 
order: 

1.  Flag  Section  (not  always  present) 

This  section  of  the  IGES  file  signals  the  format  used  to  write  the  file.  The 
absence  of  this  section  is  interpreted  to  mean  that  the  normal  ASCII  format 
was  used  when  creating  the  file. 

2.  Start  Section 

The  Start  Section  is  intended  to  be  a  human-readable  prologue  containing  com¬ 
ments  about  the  IGES  file. 

3.  Global  Section 

The  Global  Section  of  the  IGES  file  contains  information  about  the  CAD  system 
that  created  the  file  and  information  that  should  be  considered  by  the  CAD 
system  interpreting  the  file  before  the  IGES  file  is  processed. 

4.  Directory  Entry  Section 

The  Directory  Entry  Section  consists  of  one  two-line  entry  for  each  entity  de¬ 
scribed  by  the  IGES  file. 

4For  a  complete  description  of  IGES,  see  reference  7. 
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5.  Parameter  Data  Section 

The  Parameter  Data  Section  of  the  IGES  file  contains  the  geometric  information 
that  will  be  used  to  reconstruct  each  entity. 

6.  Terminate  Section 

The  goal  of  an  IGES  translation  is  to  preserve  the  geometry  and  functionality  of 
entities  in  a  CAD  drawing  or  an  IGES  file.  This  process  does  have  limits.  As  with  all 
translation,  concepts  that  can  be  expressed  precisely  in  one  language  may  not  have 
exact  equivalents  in  another  language;  conversely,  concepts  common  to  two  languages 
may  be  expressed  differently  by  each.  The  situation  applies  as  much  to  CAD  data 
as  to  natural  languages.  When  entities  have  no  direct  correspondence  between  IGES 
and  the  entities  in  a  particular  CAD  program,  the  translator  maps  them  to  similar 
constructs  that  attempt  to  preserve  as  much  data  as  possible. 

For  example,  IGES  has  no  direct  equivalent  to  an  AutoCAD  tapered  polyline.  To 
translate  a  tapered  polyline,  AutoCAD  creates  an  IGES  Composite  Curve  segment 
whose  width  is  the  average  of  the  polyline’s  starting  and  ending  widths.  Similarly, 
AutoCAD  has  no  counterpart  to  the  IGES  Parametric  Spline  Surface,  so  AutoCAD 
approximates  these  entities  with  3D  meshes. 

Many  drawings  can  be  translated  with  little  or  no  loss  of  data,  but  even  in  this 
case,  the  entities  used  to  represent  the  drawing  may  change  in  translation.  This  means 
that  IGES  is  not  fully  symmetrical:  reading  an  IGES  file  with  the  same  program  that 
was  used  to  create  it  does  not  necessarily  lead  to  a  drawing  that  is  identical  to  the 
original.  The  more  complex  the  drawing,  the  more  likely  that  information  will  have 
to  be  approximated  (this  applies  especially  to  drawings  that  are  heavily  annotated  or 
hierarchically  organized,  or  that  use  complex  three-dimensional  entities).  For  a  one¬ 
time  translation  to  or  from  IGES,  this  may  not  pose  a  great  problem.  If,  however, 
one  is  concerned  with  maintaining  drawings  that  must  be  translated  between  the  two 
systems  over  a  period  of  time,  it  is  necessary  to  be  familiar  with  the  details  of  both 
the  send  and  receive  translators  as  well  as  the  details  of  IGES  formats. 

IGES  has  been  distributed  in  successive  versions,  with  each  version  providing  ad¬ 
ditional  features  and  enhancements  to  existing  features.  Because  it  is  intended  as  a 
long-term  standard,  a  new  IGES  version  attempts  to  support  all  features  that  have 
been  officially  part  of  any  earlier  version. 

2.3.2  IGES  FOR  AUTOCAD 

The  AutoCAD  IGESIN  and  IGESOUT  commands  support  translation  of  draw¬ 
ings  to  and  from  IGES.  AutoCAD  generates  files  that  are  compatible  with  IGES  4.0. 
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It  can  successfully  read  files  that  conform  to  the  IGES  2.0,  3.0,  and  4.0  file  formats 
provided  they  employ  the  fixed-length  ASCII  form.  IGES  entities  supported  by  Au¬ 
toCAD  are  shown  in  Figures  1  and  2. 

2.3.3  IGES  FOR  ACAD 

ACAD  has  the  ability  to  transfer  a  drawing  to  other  systems  by  creating  a  file 
having  the  IGES  format.  Limitations  on  the  transfer  of  drawings  are  a  result  of 
the  limitations  of  the  IGES  translators  of  both  the  sending  and  receiving  systems. 
The  ACAD/IGES  translator  is  used  to  transfer  drawings  between  ACAD  and  other 
systems  such  as  CADAM,  CATIA  and  AutoCAD.  The  ACAD/IGES  translator  was 
created  using  the  methodology  described  [7]. 

Entities  which  are  supported  by  the  current  version  of  the  ACAD/IGES  translator 
are  listed  in  Figures  3  and  4.  Several  points  should  be  noted  when  using  the  ACAD 
system  to  translate  drawings  to  and  from  the  IGES  file  format: 

1.  When  writing  an  IGES  file  from  ACAD,  entities  which  are  confined  to  a  single 
plane  are  represented  as  an  IGES  entity  constructed  in  the  XY  plane  with 
an  associated  transformation  matrix  to  translate  the  entity  to  the  appropriate 
plane. 

2.  Spline  entities  written  from  ACAD  into  the  IGES  format  can  be  represented  as 
either  IGES  type  112  (Parametric  Spline  Curve)  or  type  126  (Rational  B-Spline 
Curve)  depending  on  the  selection  of  the  buttons  under  the  “Spline  Types:” 
section  of  the  “Write  IGES  File”  dialog  box.  If  the  button  “Parametric”  is 
selected,  all  spline  entities  will  be  written  as  type  112;  otherwise,  all  spline 
entities  will  be  written  as  type  126. 

3.  The  method  in  which  a  Surface  entity  is  written  into  IGES  format  from  the 
ACAD  system  will  depend  on  the  setting  of  the  “Surface  Type:”  buttons  in  the 
“Write  IGES  File”  dialog  box. 

4.  When  writing  ACAD  Face  entities  into  the  IGES  format,  several  IGES  entities 
are  created.  The  Face  entity  is  represented  as  type  144  (Trimmed  (Parametric) 
Surface)  which  consists  of  boundary  curves  that  are  represented  by  type  142 
(Curve  on  a  Parametric  Surface).  Boundary  curves,  in  turn,  are  created  through 
type  102  (Composite  Curve),  which  links  together  simple  entities  such  as  points, 
lines,  circles,  splines,  conics,  and  ellipses,  forming  a  closed  curve. 
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0  (Null)  Supported. 

100  (Circular  Arc)  Supported. 

102  (Composite  Supported. 

Curve) 

1 04  (Conic  Arc)  Supported:  1GESIN  approximates  this  with  a  Polyline. 

1 06  (Copious  Data)  Partially  supported:  IGESIN  discards  vectors  for  Forms 

3  and  13.  IGESOUT  translates  AutoCAD  Traces  and 
Solids  (that  are  not  extruded)  into  IGES  Simple 
Closed  Area  entities  (106  Copious  Data,  Form  63). 
This  loses  their  solid-fill  information,  but  retains 
their  area  information. 

108  (Plane)  IGESIN  uses  unbounded  planes  (Form  0)  for  clipping 

IGES  View  entities  (410).  Bounded  planes  (Forms  1 
and  -1)  are  translated  by  generating  a  Polyline  from 
the  bounding  curve. 

110  (Line)  Supported. 

1 1 2  (Parametric Spline  Partially  supported:  IGESIN  approximates  this  with  a 
Curve)  Polyline. 

114  (Parametric Spline  Partially  supported:  IGESIN  approximates  this  with  a 
Surface)  3D  Mesh. 

116  (Point)  Supported. 

118  (Ruled  Surface)  Supported. 

120  (Surface  Of  Supported. 

Revolution) 

122  (Tabulated  Supported.  See  the  section  "Extruded  Entities"  on 

Cylinder)  page  18. 


124  (Transformation 
Matrix) 

Supported. 

125  (Flash) 

Not  supported. 

126  (Rational  B-Spline 
Curve) 

Supported. 

128  (Rational  B-Spline 
Surface) 

Not  supported. 

130  (Offset  Curve) 

Not  supported. 

132  (Connect  Point) 

Not  supported. 

134  (Node) 

Not  supported. 

136  (Finite  Element) 

Not  supported. 

138  (Nodal  Displace¬ 
ment  &  Rotation) 

Not  supported. 

140  (Offset  Surface) 

Not  supported. 

142  (Curve  On  A 
Parametric 

Surface) 

Not  supported. 

144  (Trimmed 
Parametric 

Surface) 

Not  supported. 

146  (Nodal  Result) 

Not  supported. 

148  (Element  Results) 

Not  supported. 

150  (Block) 

Not  supported. 

152  (Right  Angular 
Wedge) 

Not  supported. 

154  (Right  Circular 

Not  sunnorted. 

202  (Angular 
Dimension) 

Supported. 

206  (Diameter 
Dimension) 

Supported. 

208  (Flag  Note) 

Supported. 

210  (General  Label) 

Supported. 

212  (General  Note) 

Supported. 

214  (Leader) 

Supported. 

216  (Linear 
Dimension) 

Supported. 

218  (Ordinate 

Supported. 

Dimension) 

220  (Point  Dimension)  Translated  to  AutoCAD  Block. 

222  (Radius  Supported. 

Dimension) 

228  (General  Symbol)  Supported. 

230  (Sectioned  Area)  Not  supported. 


Figure  1:  IGES  to  AutoCAD  conversion  table. 
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302  (Associativity 
Definition) 

304  (Line  Font 
Definition) 

306  (MACRO 
Definition) 

308  (Subfigure 
Definition) 

310  (Text  Font 
Definition) 

312  (Text  Display 
Template) 

314  (Color 

Definition) 

320  (Network 
Subfigure 
Definition) 

322  (Attribute  Table 
Definition) 


402  (Associativity 
Instance) 

•  404  (Drawing) 

406  (Property) 

408  (Singular 
Subfigure 
Instance) 

410  (View) 

412  (Rectangular 
Array  Subfigure 
Instance) 

414  (Circular  Array 
Subfigure 
Instance) 

416  (External 
Reference) 

418  (Nodal 

Load/Constraint) 

420  (Network 
Subfigure 
Instance) 

422  (Attribute  Table 
Instance) 

430  (Solid  Instance) 

600-699  ( 

MACRO  Instance) 

5001-9999 

(Implementor- 

Defined) 

10000-99999 

(MACRO 

Instance) 


Not  supported. 

Output  by  IGESOUT  if  necessary.  Partially  supported 
by  IGESIN. 

Not  supported. 

Supported. 

Not  supported. 

Not  supported. 

Not  supported. 

Not  supported.  * 


Not  supported.  Instead  of  creating  Attribute  Table 
Definitions,  IGESOUT  translates  AutoCAD  Attributes 
into  a  pair  consisting  of  a  General  Note  and  a  Prop¬ 
erty  entity,  and  attaches  these  as  Subfigure  references 
(see  section  2.2.4.4.2  of  the  IGES  5.1,  and  the  section 
"Attributes"  on  page  18). 

IGESIN  supports  Forms  1,  3,  4,  7,  13,  14,  15,  and  16. 
IGESOUT  supports  Form  3. 

Supported. 

Partial  support:  see  the  section  "Attributes"  on 
page  18. 

Supported. 


Supported. 

Supported,  except  that  IGESIN  does  not  support  the 
DO-DONT  flags. 


Not  supported. 


Supported.  See  the  section  "External  References"  on 
page  23  and  "External  Reference  —  416"  on  page  40. 

Not  supported. 

Not  supported. 


Not  supported. 

Not  supported. 

Not  supported. 

IGESOUT  uses  form  7901  for  Block  Attributes. 


Not  supported. 


Figure  2:  IGES  to  AutoCAD  conversion  table  (continued). 
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IGES  IGES 

Entity  Number  Description 


ACAD 

Entity 


100 

Arc 

102 

Composite  Curve 

104  Form  0-3 

Conic  Arc 

106  Form  1,2,11,12 

Copious  Data 

20,21,31-38,40,63 

108  Form  0 

Plane 

110 

Line 

112 

Parametric  Spline 

114 

Parametric  Spline  Surface 

116 

Point 

118  Form  1 

Ruled  Surface 

120 

Surface  of  Revolution 

122 

Tabulated  Cylinder 

124  Form  0 

Transformation  Matrix 

126  Form  0 

Rational  B— Spline  Curve 

128  Form  0 

Rational  B— Spline  Surface 

130 

Offset  Curve 

142 

Curve  on  a  Parametric  Surf. 

144 

(only  sub -entity  of  144) 
Trimmed  (Param.)  Surface 

202 

Angular  Dimension 

206 

Diameter  Dimension 

210 

General  Label 

212  Form  0,2-8 

General  Note 

214  Form  2,4 

Leader  (Arrow) 

216 

Linear  Dimension 

218 

Ordinate  Dimension 

220 

Point  Dimension 

222 

Radius  Dimension 

308 

Subfigure  Definition 

402  Form  1,3,4, 7, 

Associativity  Instance 

14,15 

404 

Drawing 

406  Form  15 

Name 

408 

Singular  Subfigure  Instance 

410 

View 

Arc/Circle 

any  combination  of :  point,  line,  circle 
spline,  conic,  offset 
Conic/Ellipse 

any  combination  of :  point,  line 

Plane 

Line 

Spline 

Arbitrary  Surface 
Point 

Ruled  Surface 
Surface  of  Revolution 
Ruled  Surface 
N/A 
Spline 

Arbitrary  Surface 
Spline 

any  combination  of :  point,  line,  circle, 
spline,  conic,  offset 
Face 

Angular  Dimension 
Diameter  Dimension 
Label  Dimension 
Text 

Arrow  Dimension 

Vertical/Horizontal/Parallel  Dimension 
Vertical  Call  Out/Horizontal  Call  Out/ 
Label  Dimension 

Vertical  Call  Out/Horizontal  Call  Out/ 
Label  Dimension 
Radial  Dimension 
Detail 
Views,  Sets 

Parent  View  (limit  3  per  file) 

Names  Details 
Ditto 

Auxiliary  View 


Figure  3:  IGES  to  ACAD  conversion  table. 
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ACAD 

Entity 

IGES  Entity 
Number 

IGES 

Description 

Geometric 

Point 

116 

Point 

Line 

110 

Line 

Circle/Arc 

100 

Arc  * 

Spline 

112 

Parametric  Spline  or 

Conic 

126  Form  0 

Rational  B- Spline  Curve 

104  Form  0,3 

Conic  Arc  * 

Ellipse 

104  Form  1 

Conic  Arc  * 

Plane 

108  Form  0 

Plane 

Surface  Identities: 

PolyConic,  Fillet 

114 

Parametric  Spline  Surface  or 

Ruled 

128  Form  0 

Rational  B-Spline  Surface 

118  Form  1 

Ruled  Surface  or 

114 

Parametric  Spline  Surface  or 

Surface  of  Revolution 

128  Form  8 

Rational  B-Spline  Surface 

120 

Surface  of  Revolution  or 

114 

Parametric  Spline  Surface  or 

128  Form  0 

Rational  B— Spline  Surface 

Sculptured,  Arbitrary,  Offset  114 

Parametric  Spline  Surface  or 

Face 

128  Form  0 

Rational  B~  Spline  Surface 

144 

Trimmed  (Param)  Surf,  and 

142 

Curve  on  a  Param.  Surf,  and 

102 

Composite  Curve 

Solids: 

Polyhedral 

Not  Implemented 

Volume 

Not  Implemented 

Primitive 

Not  Implemented 

Annotation: 

Text 

212  Form  0,6,7,8 

General  Note 

Dimensions: 

Horizontal,  Vertical,  Parallel  216 

Linear  Dimension  % 

Angular 

202 

Angular  Dimension  % 

Radial 

222 

Radius  Dimension  % 

Diameter 

206 

Diameter  Dimension  % 

Label,  Textline 

210 

General  Label  % 

Arrow 

214  Form  2,4 

Leader  (Arrow) 

Station  Label 

218 

Ordinate  % 

View  Call  Out 

110 

Lines  and 

212 

General  Note 

Structures: 

Detail 

308 

Subfigure  Definition  * 

Ditto 

408 

Singular  Subfigure  Instance 

Aux  Views 

410 

View 

Parent  View 

404 

Drawing 

*  written  in  the  XY  plane  with  attached  Transformation  Matrix  (124  Form  0). 

%  entity  contains  pointers  to  the  corresponding  General  Note  (212),  Witness  Lines  (106 
Form  40),  Arrows  (214),  and  Transformation  Matrix  (124  Form  0). 


Figure  4:  IGES  to  ACAD  conversion  table  (continued). 
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2.4  ACAD  GENERIC  FACET  FILE 


The  ACAD  Generic  Facet  file  is  an  ASCII  file  containing  faceted  (polygonal)  data 
as  well  as  edge,  vertex,  and  surface  normal  data.  The  files  are  named  with  the  “.facet” 
suffix.  The  convention  for  surface  normals  is  to  point  away  from  the  interior  of  the 
part.  The  triangular  facet  shape  is  the  only  one  currently  supported.  The  Generic 
Facet  format  version  V3.0  is  described  below: 

Line  A:  Revision  Date/ Time  Machine 

Revision  =  File  format  version; 

Date/Time  =  Date  and  time  of  creation; 

Machine  =  Hardware  platform  of  origin. 

Line  B:  NP 

NP  =  Number  of  parts  in  file. 

Line  C:  Part  Name 

Part  Name  =  Name  of  current  part. 

LineD:  MIRROR  (A  B  C  D) 

MIRROR  =  Mirrored  about  a  plane? 

(if  not  mirrored  about  a  plane  then  MIRROR  =  0  and  A,  B,  C,  D,  are  not 
present) 

LineE:  NV 

NV  =  Number  of  vertices  in  current  part. 

Line  F:  X  Y  Z 

(X  Y  Z)  =  3-D  Cartesian  point: 

(there  will  be  NV  copies  of  this  line). 

Line  G:  NSP 

NSP  =  Number  of  subparts  in  current  part. 
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Line  H:  Sub-Part  Name 
Sub-Part  =  Name  of  the  current  subpart. 

Line  J:  ET  NSE  NSV  EM2  VP  VN  EC 
ET  =  Type  of  element  (Triangle  =  3); 

NSE  =  Number  of  elements  in  current  subpart  —  always  ^  0; 

NSV  =  Number  of  vertices  in  the  current  subpart  ^  0  if  VP  =  1  or  VN  =  1; 

EM2  =  Parameter  set  to  1  if  2-sided  materiel  fields  have  been  defined,  and  0 
otherwise: 

When  set  to  1,  Line  N  will  have  3  material  fields,  instead  of  1.  In  these 
fields,  M  will  describe  the  material  associated  with  the  element  only,  Ml 
will  describe  the  material  associated  with  the  plus  normal  side  of  the  el¬ 
ement  and  M2  the  material  associated  with  the  minus  normal  side  of  the 
elements. 

VP  =  Parameter  set  to  1  if  vertex  parameters  are  present,  and  0  otherwise: 

Line  K  is  present  in  the  file  only  when  this  parameter  is  set  to  1. 

VN  =  Parameter  set  to  1  if  vertex  normals  are  present,  and  0  otherwise: 

Line  L  is  present  in  the  file  only  when  this  parameter  is  set  to  1. 

EC  =  Parameter  set  to  1  if  element  curvature  lines  are  present,  and  0  otherwise: 
Line  M  is  present  in  the  file  only  when  this  parameter  is  set  to  1. 

Line  K:  U  V  VID 

(U  V)  =  Parametric  vertex  coordinate  —  there  will  be  NSV  copies  of  this  line 
if  VP  =  1; 

VID  =  Vertex  ID  referencing  the  part  vertex  list. 

Line  L:  Nx  Ny  Nz  VID 

(Nx  Ny  Nz)  =  3-D  unit  vector  pointing  away  from  the  interior  of  the  part: 

(there  will  be  NSV  Copies  of  this  line  if  VN  =  1) 

VID  =  Vertex  ID  referencing  the  part  vertex  list. 
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Line  M:  Min  Max  MnVx  MnVy  MnVz  MxVx  MxVy  MxVz 

Min/ Max  =  The  principal  curvatures  computed  at  the  facet  center 

(positive  values  indicate  surface  bending  towards  the  normal) 

MnVx  MnVy  MnVz  =  A  3-D  unit  vector  pointing  in  the  direction  where  cur¬ 
vature  is  a  minimum; 

MxVx  MxVy  MxVz  =  Points  in  the  direction  of  maximum  curvature. 

(MnV  x  MxV  =  N) 

(there  will  be  NSE  copies  of  this  line  if  EC  =  1) 

Line  N:  VI  V2  M 

VI,  V2  =  Indices  of  subpart  vertices  if  NSV  ^  0,  otherwise  indices  of  part 
vertices; 

M  =  Fields  to  describe  properties  associated  with  this  element 
(there  will  be  NE  copies  of  this  line  if  ET  =  2) 

Line  N:  VI  V2  V3  M  (Ml  M2) 

VI,  V2,  V3  =  Indices  of  subpart  vertices  if  NSV  ^  0,  otherwise  indices  of  part 
vertices; 

M  (Ml  M2)  =  Fields  to  describe  material  properties  associated  with  this  element 
(there  will  be  NE  copies  of  this  line  if  ET  =  3) 

Line  N:  VI  V2  V3  V4  M  (Ml  M2) 

VI,  V2,  V3,  V4  =  Indices  of  subpart  vertices  if  NSV  ^  0,  otherwise  indices  of 
part  vertices; 

M  (Ml  M2)  =  Fields  to  describe  material  properties  associated  with  this  element 
(there  will  be  NE  copies  of  this  line  if  ET  =  4) 

Line  N:  VI  V2  V3  V4  V5  V6  M  (Ml  M2) 

VI,  V2,  V3,  V4,  V5,  V6  =  Indices  of  subpart  vertices  if  NSV  ^  0,  otherwise 
indices  of  part  vertices; 
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M  (Ml  M2)  =  Fields  to  describe  material  properties  associated  with  this  element 
(there  will  be  NE  copies  of  this  line  if  ET  =  6) 

Sections  H  thru  N  are  repeated  for  each  subpart. 

Sections  C  thru  N  are  repeated  for  each  part. 

2.5  SHIP  CAD  FILES 

NAVSEA  and  its  contractors  have  made  extensive  use  of  three-dimensional  ship 
drawings.  Engineers  have  their  personal  preferences  when  it  comes  to  CAD  soft¬ 
ware,  and  therefore  the  IGES  file  format  is  extremely  important  in  the  exchange  of 
databases  between  platforms  and  programs.  AutoCAD  binary  files  for  the  ships  listed 
in  Table  2  were  provided  to  NPS  by  NAVSEA.  The  files  were  imported  to  AutoCAD 
residing  on  NPS  platforms,  and  then  translated  to  IGES  using  the  IGES  OUT  com¬ 
mand.  Next,  the  output  file  is  read  into  ACAD,  and  an  ACAD  binary  file  written  for 
future  structure  modification.  Also,  an  ACAD  facet  file  can  be  generated  for  trans¬ 
lation  to  PATCH  input  format. 

Three-dimensional  ship  drawings  generated  by  AutoCAD  are  shown  in  Figures  5 
through  7  for  a  DDG-51,  LHA,  and  LHD,  respectively.  Figure  8  illustrates  the  level 
of  detail  that  is  typical  of  the  ship  databases.  When  performing  EM  computations 
at  HF  this  level  of  detail  is  not  required;  only  structures  on  the  order  of  0.1A  need  to 
be  represented.  (At  the  high  end  of  the  HF  band  the  frequency  is  30  MHz  and  the 
wavelength  10  meters.)  Furthermore,  meshing  the  surfaces  of  small  structures  intro¬ 
duces  a  large  additional  number  of  edges  (i.e.,  unknowns)  that  do  not  improve  the 
quality  of  the  solution,  but  dramatically  increases  the  computer  run  time.  Therefore, 
it  is  necessary  to  filter  the  geometry  model  to  remove  unnecessary  detail.  This  step 
can  be  done  either  in  AutoCAD  (before  the  IGES  file  is  written),  or  after  the  model 
has  been  imported  to  ACAD. 

3.0  PATCH 

3.1  PATCH  GENERAL  DESCRIPTION 

PATCH  [1]  is  a  FORTRAN  computer  code  that  computes  electromagnetic  scat¬ 
tering  and  radiation  based  on  a  method  of  moments  (MM)  solution  of  the  E-field 
integral  equation  (EFIE).  The  method  of  moments  reduces  the  EFIE  to  a  set  of  lin¬ 
ear  equations  that  can  be  solved  using  standard  matrix  methods.  The  number  of 
unknowns,  and  hence  the  size  of  the  matrix  equation  that  must  be  solved,  depends 
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Figure  5:  AutoCAD  generated  model  of  a  DDG-51  (without  weapons  systems  and 
antennas). 
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Figure  6:  AutoCAD  generated  model  of  a  LHA. 
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Table  2:  Ship  3D  CAD  Files 


Ship 

Imported 
in  ACAD 

Meshed 
in  ACAD 

DDG51 

X 

LHA 

X 

LHD* 

X 

DD963 

X 

X 

*  No  hull  file 


on  the  number  of  triangular  patches  that  are  used  to  represent  the  scattering  body. 

The  method  of  moments  is  rigorous;  that  is,  in  the  limit  as  the  triangles  become 
smaller,  the  method  of  moments  solution  converges  to  the  correct  value.  Unlike  a 
wire  simulation  (wire  grid  model),  the  area  between  edges  which  form  the  triangle 
facets  are  solid  material,  not  air  gaps.  Therefore,  current  truely  flows  on  the  surface 
of  the  object,  not  just  along  the  edges  of  the  triangles. 

For  an  object  with  N  edges,  PATCH  computes  a  vector  of  complex  coefficients 
7m,  m  =  1,2, ...,  N,  such  that  the  current  crossing  edge  m  is 

Jrn 

The  unit  vector  nm  is  in  a  direction  normal  to  edge  m  and  lies  on  the  surface.  Once 
the  current  coefficients  have  been  determined  using  the  method  of  moments  proce¬ 
dure,  it  is  possible  to  compute  radiation  patterns  and  scattered  fields.  A  summary  of 
the  capabilities  of  PATCH  is  given  in  Table  3.  The  details  of  modeling  the  deckedge 
antennas  are  discussed  in  [8]. 

The  application  of  MM  is  usually  limited  by  the  size  of  the  computer  available. 
Bodies  comprised  of  large  numbers  of  triangles  yield  matrices  too  large  for  the  com- 
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Table  3:  Summary  of  PATCH  Capabilities 


Arbitrary  Shape 

Open/Closed  objects 
Modelled  by  triangular  “patches” 

Variable  patch  density 
Front  end  for  graphical  composition 
Arbitrary  edge  multiplicity 
Non-orient  able  surfaces  (e.g.,  Moebius  strip)  OK 
Symmetry  planes  may  be  included 
Multiple  bodies  okay 

Surface 

Basis  functions  yield  surface  currents 
Type:  Wilton-Rau 
Free  of  line  and  point  charges 
Equivalent  Thevenin  circuits  can  be  calculated 
Lumped  and  surface  impedance  loading  possible 


Excitation 

Voltage  sources  (e.g.,  for  antennas) 

Plane  waves 

Both 

Calculated  Quantities 
Surface  currents 
Far  field  patterns 
Radar  cross  sections 

Field  calculations  at  general  observation 
points  (including  near  field  points) 

Frequency  Domain 
Pattern  loops 
Frequency  loops 


puter,  or  run  times  too  long  to  be  of  practical  use.  A  general  rule  of  thumb  for 
convergence  of  the  far  field  is  that  triangle  edge  lengths  should  not  exceed  0.1  A, 
where  A  is  the  wavelength.  The  number  of  edges  is  generally  a  close  estimate  of  the 
number  of  unknowns  that  must  be  determined.  Differences  are  due  to  the  fact  that 
some  edges  may  be  shared  by  more  than  two  facets,  and  therefore  the  number  of 
current  coefficients  associated  with  that  edge  is  more  than  two.  (This  is  referred  to 
in  the  PATCH  manual  as  the  multiplicity  of  the  edge.)  A  SGI  Indigo  II  workstation 
with  128  MBytes  of  memory  can  handle  structures  with  approximately  6000  edges. 

3.2  PATCH  INPUT  FILE  FORMAT 

The  PATCH  input  file  is  an  ASCII  file  that  contains  all  of  the  geometry  informa¬ 
tion  and  calculation  parameters.  The  file  can  be  generated  using  the  preprocessing 
code  BUILDN5  which  is  distributed  along  with  PATCH.  BUILDN5  is  capable  of  gen¬ 
erating  basic  geometrical  shapes  and  combining  them  to  yield  more  complex  shapes. 
In  addition  to  geometry  information,  BUILDN5  also  prompts  the  user  for  calculation 
information  such  as  frequency,  observation  angles,  and  excitation  conditions. 

The  user  input  data  is  appended  to  the  geometry  file  and  written  in  a  format 
that  is  recognizable  to  PATCH.  Upon  execution,  the  NPS  version  of  PATCH  (which 
has  been  modified  from  the  original)  looks  for  a  file  named  “inpatch”  in  the  current 
directory.  It  performs  the  required  calculations  and  generates  an  ASCII  output  file 
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named  outpatch  which  contains  all  of  the  input  data  in  readable  form.  An  output 
file  of  current  coefficients  is  also  generated. 

It  is  not  necessary  that  BUILDN5  be  used  to  generate  the  input  file.  If  the  user 
is  familiar  with  the  required  file  format,  “inpatch”  can  be  generated  using  any  text 
editor.  The  input  data  format  follows:5 


Line  A:  TITLE 

TITLE  =  80  character  ASCII  string. 

LineB:  NNODES,NEDGES 
NNODES  =  Number  of  nodes  in  file. 

NEDGES  =  Number  of  edges  in  file. 

Line  C:  NODE,  X(NODE),  Y(NODE),  Z(NODE) 

NODE  =  Node  index. 

X  =  X  coordinate  of  node  number  NODE. 

Y  =  Y  coordinate  of  node  number  NODE. 

Z  =  Z  coordinate  of  node  number  NODE. 

(there  will  be  NNODES  copies  of  this  line) 

Line  D:  NEDGE,  NODEl(NEDGE),  NODE2(NEDGE) 

NEDGE  =  Edge  number. 

NODE1  =  Node  number  of  first  end. 

NODE2  =  Node  number  of  second  end. 

(there  will  be  NEDGES  copies  of  this  line;  the  order  of  endpoints  is  not 
important) 

Line  E:  IGNDP(l),  IGNDP(2),  IGNDP(3) 

5This  input  sequence  is  typical  for  the  calculation  of  antenna  patterns  or  received  signals  due  to 
plane  wave  incidence.  It  does  not  necessarily  cover  all  possible  input  sequences.  See  reference  1  for 
a  complete  discussion  of  the  input  format. 
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IGNDP  =  Symmetry  planes  at  x  =  0,  y  =  0  or  z  =  0? 

0  =  no  ground  plane;  1=  infinite  PMC;  -1  =  infinite  PEC 

Line  F:  NEXCIT 

NEXCIT  =  Number  of  voltage  excitations. 

Line  G:  ITYPE 

ITYPE  =  Type  of  excitation. 

p  =  plane  wave;  v=  voltage;  b  =  both 
(there  are  NEXCIT  copies  of  this  line) 

Line  H:  If  ITYPE  =  p:  THETA,  PHI,  ETRE,  ETIM,  EPRE,  EPIM 
THETA,  PHI  =  ( 6 ,  <j>)  angle  of  incidence  in  polar  coordinates. 

ETRE, ETIM  =  Real(Ee),  Imag(^). 

EPRE, EPIM  =  Real(E^),  Imag(E^). 

Line  I:  IMAG(l),  IMAG(2),  IMAG(3) 

IMAG  =  Image  the  plane  wave  about  x  =  0,  y  =  0,  or  z  =  0? 

0  =  no  image;  1=  image  for  PMC;  -1  =  image  for  PEC 

Line  H:  If  ITYPE  ^  p:  NVOLT 
NVOLT  =  Number  of  voltage  sources. 

Line  I:  IEDGV,  IPOS,  VREAL,  VIMAG 
IEDGEV  =  Edge  number  for  face  on  which  the  voltage  source  resides. 
IPOS  =  Node  opposite  IEDGEV  for  positive  voltage  sense. 

VREAL  =  Real  part  of  impressed  voltage. 

VIMAG  =  Imaginary  part  of  impressed  voltage. 

(there  will  be  NVOLT  copies  of  this  line) 

Line  J:  NFZS 
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NFZS  =  Number  of  faces  with  nonzero  surface  impedance. 

Line  K:  IFZS,  ZSRE,  ZSIM 
IFZS  =  Face  number. 

ZSRE  =  Real  part  of  surface  impedance. 

ZSIM  =  Imaginary  part  of  surface  impedance. 

(there  will  be  NFZS  copies  of  this  line) 

Line  L:  THEV,  MTHEV,  IETHEV 

THEV  =  Thevinin  equivalent  circuit? 

If  THEV  =  .false,  then  MTHEV  =  0  and  IETHEV  =  0 
If  THEV  =  .true,  then  MTHEV  =  basis  function  index 
If  THEV  =  .true,  then  IETHEV  =  edge  index 

Line  M:  IPATT 
IPATT  =  0,  no  pattern  calculation. 

IPATT  =  1,  pattern  calculation  with  3-point  integration. 
IPATT  =  2,  pattern  calculation  with  1-point  integration. 

Line  N:  If  IPATT  ^  0:  TH1,  TH2,  NTH,  PHI,  PH2,  NPH 
TH1,  TH2  =  9  pattern  limits 
PHI,  PH2  =  <j>  pattern  limits 
NTH,  NPH  =  number  of  pattern  points  in  9  and  <f> 

Line  O:  NNFLD 

IPATT  =  Number  of  field  observation  points. 

Line  P:  If  NNFLD  ^  0:  DX,  DY,  DZ 
DX,  DY,  DZ  =  Finite  difference  increments  in  a;,  y,  z. 

Line  Q:  NNFLD  /  0:  RFLD(1,J),  RFLD(2,J),  RFLD(3,J) 
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RFLD(1-3,J)  =  z,y,  z  coordinates  of  field  point  j 
(there  will  be  NNFLD  copies  of  this  line) 

Line  R:  PRINTC 
IPATT  =  .true.,  print  a  current  table. 

Line  S:  NEDGO 

NEDGO  =  number  of  edges  in  current  table. 

Line  T:  If  NEDGO  ±  0  NEGCUR 
NEGCUR  =  Edges  in  current  table. 

(there  will  be  NEDGO  copies  of  this  line) 

Line  U:  FREQ  or  -1 
FREQ  =  Frequency. 

(terminates  when  FREQ  =  -1) 

3.3  PATCH  CODE  MODIFICATIONS 

The  PATCH  source  code  received  from  Sandia  Labs  has  been  modified  to  provide 
additional  capabilities  that  were  not  available  in  the  original  version.  Thus  PATCH 
actually  refers  to  a  collection  of  codes.  The  particular  versions  of  interest  for  this 
application  are: 

1.  PATCH2V 

There  are  no  significant  changes  between  this  version  and  the  orignal  one 
provided  by  Sandia.  The  major  change  is  the  addition  of  “facet  checking” 
as  described  below. 

2.  PATCHDF 

This  version  has  been  modified  to  specifically  compute  the  current  induced 
by  incident  plane  waves  at  ship  deckedge  antenna  locations.  The  edge 
indices  corresponding  to  the  deckedge  antenna  locations  must  be  provided. 
The  induced  currents  for  all  specified  incidence  angles  are  written  to  a 
ASCII  file  for  use  by  the  programs  RECAL  and  RMSDF.  This  version  also 
does  “facet  checking.” 
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Figure  9:  An  example  of  an  ambiguous  edge-defined  surface. 

3.3.1  INPUT/OUTPUT 

All  versions  of  PATCH  run  at  NPS  generate  a  set  of  MATLAB  “.m”  files.  The 
files  contain  the  computed  field  quantities,  and  can  be  loaded  into  MATLAB  for  plot¬ 
ting.  An  ASCII  file  of  the  computed  current  coefficients  named  “currents”  is  also 
generated.  This  is  a  duplicate  of  the  fist  that  would  occur  in  “outpatch”  if  the  user 
requested  a  listing  of  the  currents. 

3.3.2  FACET  CHECKING  SUBROUTINE 

PATCH  defines  geometrical  shapes  on  the  basis  of  edge  connections;  i.e.,  it  searches 
for  three  connected  edges  and  considers  the  enclosed  edges  to  define  a  unique  facet. 
Figure  9  illustrates  a  situation  that  occurs  frequently  in  the  meshing  of  surfaces  that 
PATCH  misinterprets.  PATCH  finds  four  triangles  (ABD,  ACD,  CDB,  and  ABC) 
when  in  fact  there  are  only  three  (ABD,  ACD,  and  CDB).  Therefore,  before  declar- 
ing  that  a  face  has  been  found,  the  new  face  should  be  checked  against  all  previously 
defined  faces  to  see  if  they  have  any  common  area. 

An  efficient  test  based  on  a  comparison  of  circles  inscribed  inside  of  the  two  tri¬ 
angles  under  consideration  is  illustrated  in  Figure  10.  If  the  inscribed  circles  overlap 
then  the  two  facets  share  area  and  the  larger  of  the  two  triangles  is  not  a  valid  face. 
This  test  can  fail  in  the  case  of  extreme  aspect  ratios  as  shown  in  Figure  11.  However, 
the  “Shell  Mesh”  parameters  can  be  set  to  avoid  this  situation. 
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Figure  10:  A  test  for  common  area  based  on  a  comparison  of  inscribed  circles. 


Two  triangles  need  to  be  tested  only  if: 

1.  they  share  a  common  edge,  and 

2.  they  lie  in  the  same  plane6. 

Referring  to  Figure  10,  the  triangle  nodes  are  Pij  where  the  first  subscript  denotes  the 
triangle  number  ( i  =  1,2)  and  the  second  the  node  number  (j  =  1,2,3).  Similarly 
the  edges  are  defined  by  Lij,  which  in  vector  form  are 

U i  =  (z;2  -  z;i)z  +  ( yi2  -  yn)y  +  (zi2  -  za)z 

Li2  -  (Xi3  -  xi2)x  +  ( yi3  -  yi2)y  +  (zi3  -  zi2)z 

Liz  =  (x{1  -  xi3)x  +  (yix  -  yi3)y  +  (zix  -  zi3)z 

Position  vectors  to  the  nodes  are: 

P ( Pij )  XijX  T  yijy  T  z^j z 

The  perimeter  of  triangle  i  is 

Ci = x:  i4i 

3=1 

6 Unless  they  form  the  open  end  of  a  pyramid  or  corner.  This  condition  will  not  be  encountered 
in  general. 
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Figure  11:  An  example  of  a  case  in  which  the  test  fails. 

In  terms  of  the  node  position  vectors  the  area  of  triangle  i  is  given  by 

Ai  =  | R(Ph)  x  R(Pi2 )  +  £(^2)  x  R(Pi3)  +  R(Pi3)  x  £(*>•,  )| 

A  normal  vector  is  the  cross  product  of  any  two  edges.  For  instance, 

Ni  —  L{\  x  Aj'2 

Finally,  the  position  vector  to  the  center  of  the  inscribed  circle  for  triangle  i  is 

Oi  =  [Li2\R(Pa)\  +  Li3\R{Pi2)\  +  4|4^'3)|]  /Ci 

The  radius  of  the  circle  is  r,-  =  Aj/C',-.  Two  circles  i  and  j  overlap  if  the  distance 
between  their  centers  is  less  than  the  sum  of  their  radii 

I  Oi  0j\  <  r,-  +  r^- 

Using  the  above  equations,  the  following  test  can  be  applied  to  determine  whether 
or  not  a  triangle  is  a  valid  face: 

1.  Loop  through  all  pairs  of  triangles.  For  two  triangles  that  share  and  edge, 

2.  see  if  they  lie  in  the  same  plane  (to  within  some  tolerance). 

3.  If  they  do,  find  out  which  triangle  has  the  smallest  area  (smallest  inscribed 
circle)  and  find  the  location  of  the  center. 
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4.  Determine  if  the  circles  overlap  by  finding  the  distance  between  center  and 
comparing  it  to  the  sum  of  the  radii. 

5.  If  the  circles  overlap  the  large  triangle  is  not  a  valid  face. 

The  facet  checking  algorithm  has  been  incorporated  into  all  of  the  PATCH  related 
codes  that  generate  a  face  list  from  the  edges.  They  include: 

All  versions  of  PATCH  (new  subroutine  GEOM  and  additional  subroutines  AXB 
and  FACETCK) 

PATCH-to-ACAD  translator,  PTA 
Input  file  checking  program,  KNIT 
MATLAB  plot  file  generating  program,  BLDMAT 
The  required  changes  to  the  PATCH  codes  to  implement  facet  checking  are: 
new  subroutine  GEOM 

additional  subroutine  AXB  (vector  cross  product) 
additional  subroutine  FACETCK 

These  codes  are  listed  in  Appendix  A.  Note  that  the  argument  list  in  GEOM  is  not 
the  same  as  that  for  the  original  version. 

4.0  VIEWING  GEOMETRY  FILES 

4.1  INTRODUCTION 

Geometry  files  can  be  viewed  using  several  different  methods.  The  information 
that  can  be  displayed  differs  in  each  case.  The  options  available  are: 

1.  DISSPLA  graphics  via  the  program  BUILDN5  (Facet,  node,  and  edge  numbers 
can  all  be  displayed.) 

2.  ACAD  (Only  facet  numbers  are  displayed.  Edge  and  node  numbers  can  be 
found  using  the  edge  connection  fist  generated  by  the  program  BLDMAT.) 

3.  MATLAB  (Only  node  and  edge  numbers  can  be  displayed.) 
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The  procedures  required  to  use  each  of  the  above  three  methods  are  described  below. 

4.2  VIEW  USING  BUILDN5 

A  PATCH  input  file  can  be  viewed  on  the  screen  or  printed  using  a  postscript 
file  by  running  BUILDN5.  The  “disp”  option  is  chosen  after  the  data  file  has  been 
read.  The  program  prompts  for  the  quantities  to  be  displayed  (nodes,  edges,  or  faces), 
the  limits  of  the  viewing  box,  and  viewing  angle.  The  DISSPLA  software  package  is 
required  for'  this  method. 

4.3  VIEW  USING  ACAD 

Geometries  can  be  viewed  directly  in  ACAD  after  they  have  been  meshed.  How¬ 
ever,  if  facet  numbers  in  “inpatch”  are  to  be  the  same  as  those  viewed  in  ACAD,  it  is 
necessary  to  convert  the  PATCH  input  file  back  to  the  “.facet”  form  using  PATCH’S 
index  ordering.  This  is  achieved  using  the  PATCH-to-ACAD  translator  PTA.  PTA 
creates  a  the  file  4  out. facet”  from  an  input  file  “in.patch”  in  which  each  face  de¬ 
fined  by  PATCH  is  written  as  a  separate  part.  Therefore,  facet  numbers  assigned  by 
PATCH  will  be  the  same  as  the  facet  numbers  viewed  in  ACAD  using  the  “Verify 
Entity”  command. 

It  is  not  possible  to  find  edge  numbers  visually  in  ACAD  with  the  current  trans¬ 
lator  software.  Edge  and  node  numbers  can  be  found  indirectly  by  noting  the  two 
face  numbers  that  are  attached  to  the  edge.  The  edge  number  can  be  determined  by 
finding  the  common  edge  in  the  connection  list  generated  by  PATCH  or  BLDMAT. 

4.4  VIEW  USING  MATLAB 

The  geometry  can  be  displayed  using  MATLAB.  First  the  geometry  file  (in  PATCH 
format)  is  converted  to  a  set  of  “.m”  files  using  the  FORTRAN  program  BLDMAT. 
After  BLDMAT  has  been  executed,  the  MATLAB  script  PLTPATCH  can  be  run. 
PLTPATCH  has  several  flags  that  control  whether  a  wire  grid,  surface,  or  surface 
with  hidden  lines  is  displayed.  Flags  can  also  be  set  to  display  edge  and  node  num¬ 
bers.  The  scale  and  view  can  be  changed  using  the  standard  MATLAB  commands. 

5.0  TRANSLATORS  AND  COMPUTER  CODES 

5.1  INTRODUCTION 

Throughout  the  course  of  this  research  several  new  computer  codes  were  written 
to  perform  data  translation  and  manipulation.  The  bulk  of  the  translator  codes  are 
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simply  subroutines  that  have  been  extracted  from  PATCH,  with  some  minor  mod¬ 
ifications.  The  functions  and  relationships  between  the  various  computer  codes  are 
illustrated  in  Figure  12.  (Names  with  a  “.x”  extension  refer  to  executable  files  rather 
than  the  FORTRAN  source  codes  which  carry  a  “.f”  extension.) 

Among  the  codes  developed  are  the  following  translators: 

1.  ACAD  to  PATCH 

2.  PATCH  to  ACAD 

3.  PATCH  to  NEC 

NEC  (Numerical  Electromagnetics  Code)  is  a  computational  EM  code  based  on  wire 
grid  models  [9]. 

5.2  ACAD-TO-PATCH  TRANSLATOR 

The  ACAD-to-PATCH  translator  ATP  converts  an  ACAD  “.facet”  file  (described 
in  Section  2.4)  to  a  PATCH  input  file  (described  in  Section  3.2).  A  fisting  of  program 
ATP  is  given  in  Appendix  B.  The  ACAD  file  defines  the  geometry  using  a  facet  table 
which  contains  the  node  coordinates  for  each  triangle.  ATP  defines  edges  for  each 
triangle  and  forms  an  edge  connection  fist.  It  uses  the  same  algorithm  that  is  used 
in  PATCH’S  subroutine  GEOM.  Therefore,  facet  checking  using  the  inscribed  circle 
method  is  also  incorporated  into  ATP. 

5.3  PATCH-TO-ACAD  TRANSLATOR 

The  PATCH-to-ACAD  translator  PTA  converts  a  PATCH  input  file  into  ACAD’s 
“.facet”  format.  There  are  actually  two  versions  of  this  translator:  PTA  and  PTF. 
PTA  writes  the  entire  ship  as  a  single  part.  Therefore,  when  loaded  into  ACAD,  the 
individual  triangles  cannot  be  manipulated  (i.e.,  deleted,  moved,  verified,  etc.).  PTF 
writes  an  ACAD  “.facet”  file  where  each  facet  is  a  part.  PTF  results  in  a  much  larger 
file  than  PTA.  Listings  are  given  in  Appendix  C. 

5.4  PATCH-TO-NEC  TRANSLATOR 

NEC  is  relatively  old  compared  to  the  EM  patch  codes  and  therefore  more  widely 
distributed.  A  crude  translator  was  written  so  that  a  PATCH  input  format  file  could 
be  run  on  NEC.  The  translator  maps  each  facet  edge  to  a  wire  segment.  This  can 
cause  problems  because  circular  wires  are  usually  represented  by  thin  strips  in  patch 
codes.  Thus  the  two  edges  of  a  thin  strip  result  in  two  closely  spaced  parallel  wires 
that  may  share  common  space  if  the  wire  radii  are  large  enough.  This  problem  is 
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Figure  12:  Computer  code  functions  and  relationships. 


probably  not  severe  enough  to  cause  NEC  to  abort,  but  will  affect  the  computed 
current  on  the  offending  edges. 

5.5  FILE  CHECKING  USING  KNIT 

The  program  KNIT  reads  a  PATCH  format  file  and  checks  it  for  duplicate  edges 
and  nodes.  This  can  occur  when  ACAD  creates  a  body  of  revolution.  For  example, 
a  cy Under  can  be  created  by  rotating  one  line  about  a  second  line.  This  is  analogous 
to  wrapping  a  sheet  of  paper  to  form  a  cylinder.  Where  the  two  paper  edges  meet, 
two  lines  are  created  by  ACAD.  The  presence  of  two  overlapping  edges  can  possibly 
cause  problems  when  running  PATCH  and  should  be  eliminated.  KNIT  removes  the 
duplicate  quantities  and  shifts  the  indices  of  subsequent  entries  in  the  node  and  edge 
tables  down  by  one. 

6.0  AUTOMESHING  PROCEDURE 

The  following  steps  are  used  to  generate  a  triangular  facet  model  of  a  structure 
that  has  been  created  in  AutoCAD  using  ACAD’s  “Shell  Mesh”  command  In  the 
following  discussion  it  is  assumed  that  the  reader  has  a  basic  knowledge  of  ACAD 
commands. 

Step  1:  Create  an  IGES  file  of  the  database  in  AutoCAD  using  IGESOUT. 

Step  2:  Import  the  IGES  file  into  ACAD  using  “Read  Other  Format”  under 
the  “File”  box. 

Step  3:  After  successfully  reading  the  IGES  file,  save  the  data  using  the  binary 
option  under  the  “File”  menu.  The  file  name  should  have  a  “,a9”  extension. 
This  file  will  serve  as  a  backup  in  the  event  that  ACAD  crashes  during  one  of 
the  subsequent  steps. 

Step  4:  Modify  the  file  to  suit  the  current  problem  objectives.  This  may  in¬ 
volve  discarding  entities  that  contain  too  much  detail.  Examples  are  window 
frames,  stairways,  and  internal  deckhouse  structures  such  as  shelves.  Overlap¬ 
ping  surfaces  must  also  be  eliminated. 

Step  5:  Once  the  ship  surface  has  been  uniquely  defined  by  non-overlapping 
surfaces,  create  faces  on  all  surfaces.  The  “Include  Surface  Boundaries”  option 
should  be  used. 

Step  6:  Create  a  “Sheet  Assembly”  from  the  faces  in  Step  5.  It  is  recommended 
that  a  large  minimum  edge  length  be  used  on  the  first  attempt  to  create  the 


33 


assembly.  Also,  “Cluster”  and  “Deviation”  parameter  values  of  50%  are  recom¬ 
mended.  Note  that  the  “Facegap”  parameter  will  impact  the  final  continuity  of 
the  surface;  i.e.,  whether  adjacent  faces  share  a  common  edge  or  a  gap  exists 
between  them.  The  continuity  of  the  surface  can  be  checked  using  the  “Verify 
Laminar  Edges”  option. 

Step  7:  Once  the  shell  mesh  has  been  generated  and  the  solid  volume  defined, 
turn  off  (blank)  all  entities  except  the  volume.  Save  the  volume  in  a  “.facet” 
file. 

Step  8:  Copy  the  “.facet”  file  from  Step  7  to  a  file  named  “in.facet.”  Run 
the  ACAD-to-PATCH  translator,  ATP.  A  PATCH  format  file  is  written  to 
“out.  patch.” 

Step  9:  The  program  KNIT  can  be  run  to  check  the  file  “out .patch”  for  dupli¬ 
cate  nodes  and  edges. 

Step  10:  Add  the  calculation  parameters  to  the  geometry  file.  This  can  be  done 
using  a  text  editor  (if  the  user  is  familiar  with  the  PATCH  input  file  format), 
or  by  running  BUILDN5. 

Step  11:  Copy  the  file  obtained  in  Step  10  to  one  named  “inpatch”  and  execute 
the  desired  version  of  PATCH. 

Here  are  several  useful  hints  that  should  save  time  reduce  errors: 

1.  The  body  surfaces  (usually  of  type  “Ruled”  or  “Net”)  should  be  displayed 
as  curved  mesh  using  U  and  V  grid  parameters  of  2.  This  allows  the  surfaces  to 
be  viewed  and  their  integrity  verified  without  undue  cluttering.  Before  creating 
the  faces  in  Step  5,  the  U  and  V  grid  parameters  can  be  changed  to  1.  Note 
that  in  many  cases  the  “Fit  Tol”  parameter  can  be  important. 

2.  Individual  surfaces  can  be  verified  using  the  “Verify”  command  with  the 
“Amb”  toggle  switched  on.  Triangles  that  have  only  two  flashing  edges  must  be 
deleted  and  recreated,  or  the  endpoints  flipped  on  one  of  the  edges.  (Caution: 
this  affects  the  integrity  of  any  facet  attached  to  the  edge  being  flipped.) 

3.  The  shell  meshing  (Step  6)  should  be  performed  on  moderately  sized  sub¬ 
sections  of  a  large  complex  target,  rather  than  attempting  to  mesh  the  entire 
structure. 

Figures  13  and  14  show  a  DD963  that  has  been  meshed  in  ACAD  using  the  above 
procedure.  Remeshing  the  surface  to  obtain  a  larger  or  smaller  grid  is  the  relatively 
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Figure  14:  DD963  ship  model  remeshed  with  a  small  grid  size. 
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simple  part  of  the  procedure  (which  starts  at  Step  6). 
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APPENDIX  A:  FACET  CHECKING  CODE 


The  following  subroutines  perform  the  facet  checking  algorithm  described  in  Sec¬ 
tion  3.3.2 


subrout  ine  geom ( datnod ,  nconn ,  nedges ,  itrak ,  nbound ,  mxf ace , 

$  nfaces,mxbdnd,nunknb) 

c  MODIFIED  TO  DO  FACET  CHECKING  —  NOTE  ARRAY  datnod  IS  REQUIRED 
c  this  subroutine  fills  nbound  with  the  faces  formed  by  the  edges  in 
c  nconn.  it  also  fills  in  the  multiplicity  factor  of  the  edge  in 
c  nconn(3,edge) .  it  returns  the  number  of  faces(nfaces) ,  and  the  number 
c  of  body  unknowns (nunknb)  which  is  equal  to  the  summation  of  the 
c  multiplicity  factors  of  the  edges  before  symmetry  planes  are  considered, 
c  itrak  is  a  work  array. 

integer  nconn (3, nedges) , nbound (3, mxf ace) , itrak (nedges) ,np(6) 

dimension  datnod(3,mxbdnd) 

nfaces=0 

c  find  faces  and  list  them  in  nbound. 
do  100  iedge=l,nedges-2 
ntrak=0 

c  look  for  all  edges  that  attach  to  edge  iedge  and  put  them  in  itrak. 
do  200  jedge=iedge+l, nedges 
do  20  i-1,2 
do  21  j=l ,2 

if (nconn(i, iedge) . eq.nconn(j , j edge)) then 
c  we  have  found  an  edge. 


ntrak=ntrak+l 
itrak (ntrak) = j  edge 
goto  200 

endif 

21 

continue 

20 

continue 

200 

continue 

c  find  all  pairs  of  edges  that  form  a  face  with  iedge. 
do  300  jedge=l,ntrak-l 

do  301  kedge=jedge+l,ntrak 
do  30  j=l,2 
do  31  k=l,2 

c  if  the  2  faces  in  itrak  have  a  common  point  and 
c  the  common  point  is  not  in  common  with  the  iedge... 

if ( (nconn ( j , itrak ( j  edge ) ) . eq. nconn (k, itrak (kedge))) .and. 
$(nconn(j  ,itrak(jedge))  .  ne.  nconn  (1,  iedge)  .and.nconn(j  ,  itrak  (j  edge)) 
$ .  ne .  nconn(2 ,  iedge)  )  )  then 
if (nfaces .eq.0)  then 
c  if  this  is  the  first  face  save  it 
nfaces=nfaces+l 
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c  put  the  face  into  nbound. 

nbound ( 1 , nf  ace  s ) = i edg  e 
nbound(2 ,nf aces )=itrak ( j  edge) 
nbound ( 3 , nf  ace  s )  = it rak (k  edge ) 
c  increment  the  multiplicity  factor  of  the  edges. 

nconn(3 , iedge) =nconn(3 , iedge) +1 
nconn ( 3 , it rak ( j  edge ) ) =nconn ( 3 , itrak ( j  edge ) )  + 1 
nconn  ( 3 , it rak (kedge ) ) =nconn (3 , itrak (kedge ) ) + 1 
go  to  311 
endif 

c  if  this  is  not  the  first  face,  check  to  see  if  it  overlaps  with  any 
c  previously  found  face. 
kl=iedge 
k2=itrak( jedge) 
k3=itrak (kedge) 
np(l)=nconn(l ,kl) 
np(2)=nconn(2,kl) 
np ( 3 ) =nconn ( 1 , k2) 
np(4)=nconn(2,k2) 
np(5)=nconn(l,k3) 
np(6)=nconn(2 ,k3) 
c  find  the  three  unique  points 
nodel=np(l) 
node2=np(2) 
do  603  ii=3,6 
npt=np(ii) 

if ((npt .ne.np(l)) .and. (npt .ne .np(2) ))  then 
c  must  be  the  third  node 
node3=npt 
go  to  602 
endif 

603  continue 

602  continue 

c  node  coordinates  of  the  first  face 
x 1 1 =datnod ( 1 , node 1 ) 
y 1 1 =dat nod ( 2 , node 1 ) 
zll=datnod(3,nodel) 
xl2=datnod(l ,node2) 
y 12=datnod( 2 , node2) 
z 1 2=datnod ( 3 , node2 ) 
xl3=datnod(l,node3) 
yl3=datnod(2,node3) 
zl3=datnod( 3 , node3 ) 
isum=0 

do  39  kface=l,nf aces 

c  nodes  of  the  face  number  kface  to  check  against 
il=nbound(l,kface) 
np(l)=nconn(l,il) 
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np(2)=nconn(2,il) 
i2=nbound (  2 ,  kf  ace ) 
np(3)=nconn(l,i2) 
np(4)=nconn(2,  i2) 
i3=nbound(3, kface) 
np(5)=nconn(l,i3) 
np(6)=nconn(2, i3) 
c  find  the  three  unique  points 
nodel=np(l) 
node2=np(2) 
do  613  ii=3,6 
npt=np(ii) 

if ((npt .ne.np(l)) .and. (npt .ne.np(2)))  then 
c  must  be  the  third  node 
node3=npt 
go  to  612 
endif 

613  continue 

612  continue 

c  node  coordinates  of  the  second  face 
x2 1 =dat nod ( 1 , node 1 ) 
y 2 1 =datnod ( 2 , node 1 ) 
z2 1 =dat nod ( 3 , node 1 ) 
x22=dat nod ( 1 , node2 ) 
y22=datnod(2,node2) 
z22=datnod(3,node2) 
x23=datnod(l ,node3) 
y23=datnod(2,node3) 
z23=datnod(3,node3) 
c  see  if  the  triangles  overlap 

call  facetck(xll,yll,zll,xl2,yl2,zl2,xl3,yl3, 

&  z!3 , x21 , y21 , z21 , x22 ,y22 , z22 , x23 , y23 , z23 , if lag) 

isum=isum+iflag 

c  if  iflag.eq.O  there  is  no  overlap  so  continue  checking 
if (iflag.eq.O)  go  to  38 

c  if  iflag.eq.l  do  not  include  the  face  being  tested;  no  need  to 
c  check  any  more  faces  since  this  face  is  being  discarded 
if (iflag.eq.l)  go  to  311 

c  if  iflag.eq.2  keep  the  face  being  tested  and  discard  face  number  kface 
if (if lag. eq. 2)  then 

c  reduce  the  multiplicity  factor  of  the  discarded  face  edges  by  one 
c  face  number  kface  has  edges  il,i2,i3 

nconn ( 3 , i 1 ) =nconn (3 , il ) - 1 
nconn(3 ,i2)=nconn(3,i2)-l 
nconn(3,i3)=nconn(3,i3)-l 
nbound ( 1 , kf  ace ) =iedge 
nbound ( 2 , kface ) =itr ak ( j  edge ) 
nbound ( 3 , kface ) =itrak (kedge ) 

c  increase  the  multiplicity  factor  of  the  new  face  edges  by  one. 
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nconn(3,kl)=nconn(3,kl)+l 
nconn(3  ,k2)=nconn(3,k2)+l 
nconn(3,k3)=nconn(3,k3)+l 
c  no  need  to  check  the  remaining  faces, 
go  to  311 

endif 

38  continue 

39  continue 

c  made  it  all  the  way  through  —  add  this  edge 
nfaces=nf aces+1 
c  put  the  face  into  nbound. 

nbound(l,nfaces)=iedge 
nbound ( 2 , nf  aces )  =  it r ak ( j  edge) 
nbound (3 , nf aces ) =itrak (kedge) 
c  increment  the  multiplicity  factor  of  the  edges. 

nconn(3 , iedge) =nconn(3 , iedge) +1 

nconn ( 3 , itrak ( j  edge ) ) =nconn (3 , itrak ( j  edge ) )  + 1 

nconn  ( 3 , itrak (kedge ) ) =nconn (3 , itrak (kedge ) ) + 1 

goto  311 
endif 

31  continue 

30  continue 

301  continue 

311  continue 
300  continue 
100  continue 

c  at  first  nunknb=~nedges .  this  would  be  incremented  by  3  for  each 
c  face  that  was  found,  therefore  nunknb=-nedges+3*nf aces . 
nunknb=3*nfaces-nedges 
return 
end 

subroutine  iacetck(xll,yli ,zil Jxl2,yl2,zl2,xl3Jyl3,zl3, 
ft  x21 , y21 , z2i , x22 ,y22 ,z22 , x23 , y23 , z23 , if lag) 

c  check  to  see  if  two  faces  overlap  based  on  the  common 
c  surface  area  of  circles  inscribed  in  the  two  triangles 
iflag=0 

c  only  test  with  triangles  that  lie  in  the  same  plane 
c  the  condition  is  that  the  mags  of  the  components  of  the 
c  cross  product  of  the  normals  must  be  <  eps 
eps=l. e-2 

c  notation:  first  number  in  a  variable  refers  to  face 
c  second  number  refers  to  node  or  edge 

c  x,y,z  components  of  the  three  edge  vectors  of  triangle  1 
egllx=xl2-xll 
eglly=yl2-yll 
egllz=z!2-zll 
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egl2x=xl3-xl2 
egl2y=y!3-yl2 
egl2z=zl3-zl2 
egl3x=xli-xl3 
egl3y=yli-yl3 
egl3z=zll-zl3 
c  edge  lengths  for  face  1 

egllm=sqrt (egllx**2+eglly**2+egllz**2) 
egl2m=sqrt (egl2x**2+egl2y**2+egl2z**2) 
egi3nRsqrt (egl3x**2+egl3y**2+egl3z**2) 
c  x,y,z  components  of  the  three  edge  vectors  of  triangle  2 
eg21x=x22-x21 
eg21y=y22-y21 
eg21z=z22-z21 
eg22x=x23~x22 
eg22y=y23-y22 
eg22z=z23-z22 
eg23x=x21-x23 
eg23y=y21-y23 
eg23z=z21-z23 
c  edge  lengths  for  face  2 

eg21m=sqrt(eg21x**2+eg21y**2+eg21z**2) 
eg22m=sqrt(eg22x**2+eg22y**2+eg22z**2) 
eg23m=sqrt(eg23x**2+eg23y**2+eg23z**2) 
c  unit  vectors  normal  to  each  face 

call  AxB(egllx, eglly ,egllz,egl3x,egl3y , egl3z,vnlx,vnly, vnlz) 

vmagl=sqrt(vnlx**2+vnly**2+vnlz**2) 

vnlx=vnlx/vmagl 

vnly=vnly/vmagl 

vnlz=vnlz/vmagl 

call  AxB ( eg21x , eg2 ly , eg21z , eg23x , eg23y , eg23z , vn2x , vn2y ,  vn2z ) 

vmag2=sqrt(vn2x**2+vn2y**2+vn2z**2) 

vn2x=vn2x/vmag2 

vn2y=vn2y/vmag2 

vn2z=vn2z/vmag2 

call  AxB ( vnlx , vnly ,  vnlz , vn2x , vn2y , vn2z , vnx , vny , vnz) 
c  if  I vnx I  and  I vny I  and  I vnz I  are  sufficiently  small  then  these  two 
c  faces  can  be  considered  coincident 

if ( (abs(vnx) .lt.eps) .and. Cabs (vny) .It . eps) .and. (abs(vnz) .It . 

&  eps))  then 

c  areas  of  the  two  triangles 

call  AxB(xll ,yil ,zll ,xl2,yl2,z!2,al ,bl , cl) 

call  AxB(xl2tyl2,zl2#xl3,yl3,zl3,a2,b2,c2) 

call  AxB(xl3>yl3,zl3,xll,yll,zll,a3,b3, c3) 

areal=sqrt ( (al+a2+a3) **2+ (bl+b2+b3) **2+ ( cl+c2+c3) **2) /2 . 

call  AxB(x21,y21,z21,x22,y22,z22,al,bl,cl) 

call  AxB (x22 , y22 , z22 , x23 , y 23 , z23 , a2 , b2 , c2) 

call  AxB (x23 , y23 , z23 , x2 1 , y21 , z2 1 , a3 , b3 , c3) 

area2=sqrt C(al+a2+a3)**2+(bl+b2+b3)**2+(cl+c2+c3)**2)/2. 
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c  compute  perimeters 

periml=egl Im+egl2m+egl3m 
per im2= eg2 lm+  eg22m+  eg23m 

c  find  position  vectors  to  the  centers  of  each  inscribed  circle 
clx= ( egl2m*xl l+egl3m*xl2+egl lm*xl3 ) /per iml 
cly=(egl2ra*yll+egl3m*yl2+egllm*yl3)/periml 
Clz=(egl2m*zll+egl3m*zl2+egllm*zl3)/periml 
c2x= ( eg22m*x21+eg23m*x22+eg21m*x23)/perim2 
C2y=(eg22m*y21+eg23m*y22+eg21m*y23)/perim2 
c2z= ( eg22m*z21+eg23m*z22+eg2im*z23 ) /per im2 

c  distance  between  centers 

dl2=sqrt((clx-c2x)**2+(cly-c2y)**2+(clz-c2z)**2) 

c  radii  of  circles  are  axe as /perimeters 
radl=areal/periml 
rad2=area2/perim2 

c  if  dl2  <  radl+rad2  then  there  is  overlap 
if (dl2 . It . (radl+rad2) )  then 
iflag=l 

if (rad2 . gt . radl)  iflag=2 
endif 
endif 
return 
end 

subroutine  AxB ( al , a2 , a3 , bl ,b2 , b3 , c 1 , c2 , c3) 

c  cross  product  of  two  vectors 
cl=a2*b3-a3*b2 
c2=a3*bl-al*b3 
c3=al*b2-a2*bl 
return 
end 
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APPENDIX  B:  ACAD-TO-PATCH  TRANSLATOR  CODE 

The  ACAD-to-PATCH  translator  code  ATP  is  described  in  Section  5.1.  It  con¬ 
verts  an  ACAD  “.facet”  file  to  a  PATCH  input  file. 


c  program  atp.f  (version  3  —  allows  subparts  &  uses  face  checking) 
c 

c  "acad  to  patch"  translator  ENHANCED  VERSION 
c 

c  this  program  reads  a  file  named  "in. facet"  then  reformats  the  data 
c  and  writes  it  to  the  file  named  "out. patch"  which  can  be  read  into 
c  "buildn5.f"  as  a  geometry  file,  characteristics  of  ACAD  facet  file: 
c  »  a  vertex  table  is  given  for  each  part  (e.g.,  each  face  of  a  cube 
c  is  considered  a  part) 

c  »  each  part  is  assigned  a  name 

c  »  node  index  is  reinitialized  for  each  part  and  nodes  common  to 
c  several  parts  are  present  more  than  once  (run  "knit.f"  on 

c  "out. patch"  to  eliminate  duplicate  edges) 

c  »  a  vertex  connection  list  is  given  for  each  part  (face  index  and 
c  its  three  vertex  indices) 

c  »  material  parameters  are  included 

c  »  multiple  parts  allowed,  and  each  part  can  have  several  subparts 
c 

c  x,y,z  are  node  coordinates  (index  is  facet  number) 
c  datnod  and  nconn  are  same  as  in  patch. f 

dimension  node (3, 6000) ,nbound(3,6000) ,nv(6000) , part (6000) 
dimension  datnod (3, 6000) , nconn(3, 6000) , istart (200) 
dimension  tmpdat (3 ,6000) ,icount(6000) , indsum(6000) 
dimension  ivmin(6000) ,ivtx(6000,500) ,iskip(6000) ,indx(6000) 
integer  tmpnod(3,6000) 
character*80  title 
character*20  part,subpt 

c  distances  less  than  eps  are  considered  the  same 
eps=0 . 

c  iverb=0  is  verbose  mode  —  progress  displayed 
iverb=0 

open(l,file=' in. facet  ' , status='old' ) 

xmin= 1 . e6 

ymin= 1 . e6 

zmin=l .e6 

xraax=-l .e6 

ymax=-l . e6 

zmax=-l . e6 

1  format (a80) 

2  format (a20) 
read(l,l)  title 
write(6,*)  'title: *  , title 
read(l,*)  nparts 
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c  keep  track  of  total  number  of  faces  and  vertices 
nvtl=0 
nftl=0 

do  1000  npart = 1 , npart s 
read (1,2)  part(npart) 
read(l,*)  nmir 

if (nmir.ne.O)  write (6,*)  'error:  structure  is  mirrored' 
read(l,*)  nv(npart) 
nvert s=nv (npart ) 
if (iverb.eq.O)  then 

write(6,+)  'part  number:  ', npart 
write(6,*)  '  name:  ', part (npart) 

write(6,*)  '  vertices:  ',nv(npart) 
endif 

c  read  node  table  for  this  subpart 
do  10  n=l, nvert s 
read(l,*)  xx,yy,zz 
nn=nvtl+n 
tmpdat (l,nn)=xx 
tmpdat (2 , nn) =yy 
tmpdat(3,nn)=zz 
dat nod ( 1 , nn) =tmpdat ( 1 , nn) 
dat  no  d  (  2 ,  nn  )  =  tmpdat  (  2 ,  nn  ) 
datnod  (  3 ,  nn)  =tmpdat  (  3 ,  nn) 
xmax=amaxl  (xmax ,  xx) 
ymax=amaxl  (ymax ,  yy) 
zmax=amaxl (zmax  ,zz) 
xmin=aminl  (xmin ,  xx) 
ymin=aminl  (ymin ,  yy  ) 
zmin=aminl  (zrain ,  zz) 

10  continue 
c  loop  through  subparts 
read(l,*)  nspts 
if (iverb.eq.O)  write(6,*) 

&  'number  of  subparts  in  current  part:  ', nspts 
do  15  kn=l, nspts 
read(l,2)  subpt 
if (iverb.eq.O)  then 

write(6,*)  'subpart  number:  ' ,kn 
write (6,*)  '  name:  ', subpt 

endif 

c  read:  el  type,  no.  faces,  no.  vertices,  em2,  vp,vn,ec 
c  restrictions:  em2=0  (one-sided  properties) 

c  vp=0  (vertex  parameters) 

c  vn=0  (no  normals  present) 

c  ec=0  (no  curvature  lines) 

read(l,*)  nsides ,nf aces,nsv,nem2,nvp,nvn,nec 

c  display  if  there  are  problems: 
if (nsides .ne.3) 
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&  write (6,*)  'nontriangular  facet  encountered ' 
if (nsv.ne .0)  write(6,*)  problem:  nsv  is  not  zero' 
if (nem2.ne.O)  write(6,*)  'problem:  EM2  is  not  zero' 
if (nvp.ne.O)  write(6,*)  'problem:  VP  is  not  zero' 

if (nvn.ne.O)  write(6,*)  'problem:  VN  is  not  zero' 

if (nec.ne.O)  write(6,*)  'problem:  EC  is  not  zero' 

c  if  this  is  a  part  nodes=nverts ;  if  this  is  a  subpart  nodes=nsv 

do  20  n=l,nfaces 
nn=n+nftl 

c  read  vertex  connection  list 

read(l,*)  ndel ,nde2,nde3,ndum 
tmpnod  ( 1 ,  nn)  =nde  1+nvt  1 
tmpnod(2,nn)=nde2+nvtl 
tmpnod (3 , nn) =nde3+nvtl 
node  ( 1 ,  nn) =tmpnod( 1 , nn) 
node(2,nn)=tmpnod(2,nn) 
node (3  ,nn) =tmpnod(3 ,nn) 

20  continue 

nftl=nftl+nfaces 
15  continue 

nvtl=nvt 1+nv (npart ) 

1000  continue 

nverts=nvtl 

nfaces=nftl 

write(6,*)  'total  number  of  vertices, faces  read=' ,nverts,nfaces 
c  find  duplicate  vertices  and  save  their  indices  in  array  iold. 
c  icount(i)  counts  the  number  of  times  vertex  i  occurs 
do  75  ivl=l,nverts 
xl=tmpdat ( 1 , ivl ) 
y l=tmpdat (2 , ivl ) 
zl=tmpdat(3, ivl) 
icount(ivl)=l 

c  keep  track  of  vertex  number  ivl  and  its  duplicates 
c  ivtx(node  number,  occurace  number,  duplicate  index) 
ivtx(ivl , icount (ivl) )=ivl 
do  73  iv2=i,nverts 
if (ivl.ne.iv2)  then 
x2=tmpdat ( 1 , iv2) 
y2=tmpdat (2 , iv2) 
z2=tmpdat (3 , iv2) 
dx=abs(xl-x2) 
dy=abs(yl-y2) 
dz=abs(zl-z2) 

if (Cdx.lt . eps) .and. (dy.lt.eps) .and. (dz.lt . eps))  then 
icount (ivl) =icount (ivl) +1 
ivtx (ivl , icount ( ivl) ) =iv2 
endif 
endif 

73  continue 
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75  continue 

c  for  each  vertex  find  the  smallest  index 
do  78  iv=l,nverts 
i vmin ( i v ) =n v  er t  s + 1 
do  78  ii=l,icount(iv) 

ivrain ( iv) =min ( ivmin ( iv) , ivtx( iv , ii) ) 

78  continue 

c  find  duplicate  vertices 
irem=0 

do  76  ii=l,nverts 

c  if  minimum  index  is  .It.  ii  this  is  a  duplicate 
if (ivmin(ii) .It .ii)  then 
irem=irem+l 
iskip ( irem) =ii 
endif 

76  continue 

c  order  the  indices  to  be  removed  from  lowest  to  highest 
idx=nverts 
do  85  il=l,irem 
do  83  i2=l,irem-il+l 
indx ( i 1 ) =min ( idx , iskip ( i2 ) ) 

83  continue 

idx=indx(il)+l 
85  continue 

if (iverb.eq.O)  then 
do  190  ix=l,irem 

190  write(6,*)  'remove  ix, iskip (ix) , indx(ix) 

endif 

c  set  all  node  indices  to  minimum  value 
do  95  nn=l,nfaces 
do  95  is=l,3 
iv=tmpnod ( is , nn) 
node ( is , nn) =ivmin( iv) 

95  continue 

do  96  nn=l,nfaces 
do  96  is=l,3 

tmpnod(is,nn)=node(is,nn) 

96  continue 

c  remove  duplicate  nodes  and  shift  remaining  nodes  down  one  for 
c  each  previous  node  removed 
do  79  ir=l,irem 
do  79  iv=indx(ir) ,nverts“ir 
datnod ( 1 , i v) =datnod ( 1 # iv+1 ) 
dat nod (  2 ,  i  v  ) =datnod ( 2 , iv+ 1 ) 
datnod ( 3 , i v ) =datnod ( 3 , iv+ 1 ) 

79  continue 

nv  er t  s =nv  er t  s - ir em 

c  check  each  node  of  each  face  to  see  how  many  previous  nodes 
c  have  been  removed  (to  determine  the  number  of  steps  to  decrement 
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c  the  index) 

do  90  nn=l,nfaces 
do  90  is=l,3 
idx=tmpnod(is ,nn) 

c  count  vertices  .It.  idx  that  have  been  removed 

c  nodes  greater  than  iskip(ir)  drop  one;  those  less  than  iskip(ir)  stay- 
iter  =0 

do  93  ir=l,irem 

c  if (idx.eq. indx(ir))  go  to  90 

*io  93  ii=l,idx 

if ( ii . eq . indx ( ir ) )  iter=it er+1 
93  continue 

node  (  is ,  nn)  =tmpnod  ( i  s ,  nn) -it  er 
90  continue 

write(6,*)  'final  number  of  vertices,  faces=' ,nverts,nfaces 
c  rescale  data  if  desired 

c  write(6,*)  'rescale  data?  (0=yes/l=no) ' 

c  read(5,*)  ans 

ans=l 

if(ans.eq.O)  then 
write(6,*)  ' xraax , xmin= ' ,xmax,xmin 
write(6 ,*)  ' ymax , ymin= ' ,  ymax , ymin 
write(6,*)  'zmax,zmin='  ,zmax,zmin 
write (6,*)  'enter  scale  factor' 
read(5,*)  fac 
do  50  n=l,nverts 
datnod ( 1 , n) =f  ac*datnod ( 1 , n) 
dat nod ( 2 , n) =f  ac*dat nod ( 2 ,  n) 
datnod (  3 ,  n)  =f  ac*datnod ( 3 , n) 

50  continue 
endif 

c  generate  edge  connection  list  using  "brute  force"  checking 
if (iverb. eq.0) 

&  write (6,*)  'start  to  generate  edge  connection  list' 
c  for  face  n  check  pairs  of  vertices  to  see  if  they  have  been  assigned 
c  an  edge  number,  exception  is  n=l. 
nconn (1,1) =node (1,1) 
nconn(2,l)=node(2, 1) 
nconn(l ,2)=node(2, 1) 
nconn(2 , 2) =node (3 , 1 ) 
nconn(l ,3)=node(l, 1) 
nconn(2,3)=node(3, 1) 
nedges=3 
ncount=3 
do  60  n=2,nfaces 

c  check  all  edges  to  see  if  the  three  edges  of  the  current  face  n 
c  have  already  been  assigned  an  index 
do  59  ns=l,3 
if(ns.eq.l)  then 
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npl=node(l,n) 
np2=node(2,n) 
elseif (ns. eq.2)  then 
npl=node(2,n) 
np2=node(3,n) 
else 

npl=node(3,n) 
np2=node  ( 1 ,  n) 
endif 

if (iverb.eq.O)  then 

write (6,*)  ' checking  edge  ',ns,'  of  face  ' ,n 
endif 

do  65  ne=l,nedges 

c  check  edge  ns  of  triangle  n  against  edge  ne 
mpl=nconn(l,ne) 
mp2=nconn(2 ,ne) 

if (((mpl .eq.npl) .and. (mp2 . eq.np2) ) .or. 

&  ((mp2. eq.npl) .and. (mpl. eq.np2) ))  then 

c  edge  is  a  duplicate  and  therefore  no  need  to  go  further 
go  to  66 
endif 

65  continue 

c  made  it  all  the  way  through  —  must  be  new  edge 
if (iverb.eq.O)  then 

write(6,*)  'new  edge  ',npl,np2 

write(6 , *)  'number  of  edges  that  have  been  defined  is  ' ,ncount 
endif 

nc ount =nc ount + 1 
nc onn ( 1 , nc ount ) =np 1 
nconn ( 2 ,ncount ) =np2 
nedges=ncount 

66  continue 

59  continue 

60  continue 

do  70  n=l,nedges 
70  nconn(3,n)=-l 

if (iverb.eq.O)  write(6,*)  'edge  connection  list  generated' 
c  write  reformated  data  to  file  "out .patch11 
open(2,file='out .patch') 
write(2,*)  title 
write (2,*)  nverts ,nedges 
do  151  n=l,nverts 

151  write(2,*)  n,datnod(l,n) ,datnod(2,n) ,datnod(3,n) 
do  200  n=l,nedges 

200  write(2,*)  n,nconn(l ,n) ,nconn(2 ,n) 

c  write  new  data  in  MATLAB  files 
open(12,f ile='xpts.m') 
open(13,f ile='ypts .m') 
open ( 14, f ile='zpts .m' ) 


49 


openClSjfile^endl.m*) 
open  ( 16 ,  f  ile= 9  end2 .  m 1 ) 

98  format  (i5) 

do  140  n=l,nverts 
write(12, 100)  datnod(l,n) 
write(13, 100)  datnod(2,n) 
write(14, 100)  datnod(3,n) 

140  continue 

100  format(f 15.4) 

do  150  n=l, nedges 
write(15, 101)  nconn(l,n) 
write(16,101)  nconn(2,n) 

150  continue 

101  format (i5) 

call  geom (datnod , nconn ,  nedges ,  indsum ,  nbound , mxf ace , nf aces , 

$  mxbdnd,nunknb) 
istart(l)=l 
istart(2)=nfaces+l 

call  pmtbnd  (nconn ,  nbound ,  istart ,  1 ,  nedges , nf aces , 1 ) 
c  get  body  parameters 

call  bodpar (datnod , nconn, nbound, nverts , nedges ,nf aces ,nunknb) 

998  continue 
stop 
end 

subroutine  pmtbnd  (nconn,  nbound,  istart ,  i, nedges  ,nf  aces  ,nbodys) 

c  this  subroutine  prints  the  edges  and  the  vertices  of  each  face, 
c  input: 

c  nconn  has  the  vertices  and  the  multiplicity  factor  for  each  edge, 
c  nbound  has  the  edges  for  each  face, 
c  istart  has  the  beginning  faces  for  each  body, 
c  i  is  the  present  body, 
c  nedges  is  the  total  number  of  edges, 

c  nfaces  is  the  total  number  of  faces, 

c  nbodys  is  the  total  number  of  bodys. 

integer  nconn(3,6000) ,nbound(3,6000) , istart (200) 
integer  nverts (6000, 3) 
open(3,file=,facelist O 
open(9 , f ile= 'f acedat * ) 
open(20,f  ile=,nodel  .mO 
open(21  ,f  ile=  ,node2  .m* ) 
open(22,file=,node3.mO 
c  loop  through  the  faces  of  this  body. 

do  10  il0=istart (i) , istart (i+l)-l 

call  facvtx(nconn, nedges, nbound(l , ilO) , nbound ( 2, ilO) , 

>nbound(3 , ilO) ,nvl ,nv2 ,nv3) 

write (3 , 98) ilO , nbound  ( 1 , ilO) ,nbound(2, ilO) ,nbound(3, ilO) ,nvl,nv2 
>,nv3 
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c  write  face  number  and  vertices  to  "facedat" 
write(9,*)il0,nvl,nv2,nv3 
write(20, 100)  nvl 
write (21, 100)  nv2 
write(22, 100)  nv3 
nverts(ilO, l)=nvl 
nverts(il0,2)=nv2 
nverts ( ilO , 3) =nv3 
10  continue 

98  fonnat(lx, 'face' ,iS, ’  has  edges’ ,3i5, '  with  vertices' ,3i5) 

100  format (i6) 

return 
end 

subroutine  bodpar (datnod,nconn,nbound,nnodes ,nedges ,nf aces ,  nunknb) 

c  this  subroutine  computes  the  body  parameters 
c  input: 

c  datnod(i,j)  i=l,2,3  are  the  x,y,z  coordinates  of  the  jth  node, 
c  nconn(3,  j) :  edge  j  runs  from  node  nconn(l,j)  to  node  nconn(2,j) 
c  nconn(3, j)  is  the  multiplicity  of  the  jth  edge, 

c  nbound(i,  j ) :  contains  the  ith  edge  of  the  jth  face  i=l,2,3. 

c  nnodes  =  the  number  of  body  nodes, 

c  nedges  =  the  number  of  edges, 

c  nfaces  =  the  number  of  faces, 

c  nunknb  =  the  number  of  body  unknowns, 

c  output : 

c  avedge  =  the  average  edge  length (met ers**2)  including  multiplicity, 

c  edgemx  =  the  maximum  edge  length(meters) . 

c  mxedge  =  the  edge  number  of  the  edge  with  length  edgemx. 
c  edgemn  =  the  minimum  edge  length (meters) . 
c  mnedge  =  the  edge  number  of  the  edge  with  length  edgemn. 
c  tarea  =  the  surface  area  of  the  scatter (meters**2) : for  thin 
c  structures  only  one  side  is  considered  in  the  surface  area, 

c  avarea  =  the  average  area  of  the  faces. 

c  mxarea  =  the  number  of  the  face  with  the  maximum  area(areamx) . 

c  mnarea  =  the  number  of  the  face  with  the  minimum  axea(areamn)  . 

c  ratio  =  the  minimum  height  to  base  ratio  over  all  faces, 

c  mnrtio  =  the  face  number  that  has  a  height  to  base  ratio  of  'ratio', 

dimension  datnod(3,6000) 
integer  nconn(3,6000) ,nbound(3,6000) 

common/params/ avedge,  edgemx,  mxedge ,  edg  emn,  mnedge ,  tarea,  avarea, 
$mxarea,  mnarea,  areamx ,  areamn,  ratio ,  mnrtio 
common/mchval/ valmax ,  valrain 
c  save  /params/ 

c  the  following  line  is  a  statement  function. 

size (x , y , z) =sqrt (x*x+y*y+z*z) 
c  initialization. 
sedgl=0 
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edgemx=valrain 
edg emn= valraax 
do  20  ie=l,nedges 
mult=nconn(3 , ie) 
nl=nconn(l,ie) 
n2=nconn(2, ie) 
x=datnod ( 1 , n2 ) “datnod ( 1 , nl ) 
y =dat nod ( 2 , n2 ) -datnod ( 2 , nl ) 
z=dat nod ( 3 , n2 ) -datnod ( 3 , nl ) 
edgl=size (x , y , z) 
s  edgl =s  edgl+mnlt  *  edgl 
if  (edgl .  gt .  edgemx)tlien 
edg emx= edgl 
mxedge=ie 
endif 

if  (  edgl .  It .  edgemn)  then 
edg emn= edgl 
mnedge=ie 
endif 

20  continue 

avedge= s  edgl/ nunknb 

c  compute  tarea,avarea,mnarea, areamn,mxarea,areamx, ratio ,and  mnrtio. 
ratio=valmax 
areamx=valmin 
areamn=valmax 
taxea=0. 

do  40  iface=l,nfaces 
is l=nbound( 1 , if ace) 
is2=nbound(2 , if ace) 
is3=nbound(3 , if ace) 

call  facvtx(nconn,nedges ,isl,is2,is3,nvl ,nv2,nv3) 

call  vt xcrd ( datnod , nnodes , nv 1 , nv2 , nv3 , xl , x2 , x3 , y 1 , y2 , y3 , z 1 , z2 , z3 

» 

xlrax3=xl-x3 

ylmy3=yl-y3 

zlmz3=zl-z3 

x2mx3=x2-x3 

y2my3=y2-y3 

z2mz3=z2-z3 

x2mxl=x2-xl 

y2myl=y2-yl 

z2mzl=z2-zl 

c  compute  area  of  face  by  taking  the  cross  product  of  two  edge  vectors. 
vx=ylmy3*z2raz3-zlmz3*y2my3 
vy=zlmz3*x2mx3-xlmx3*z2mz3 
vz=xlmx3*y2my3-ylmy3*x2mx3 
area=.5*size(vx,vy,vz) 

c  compute  the  square  of  the  lengths  of  each  side. 
rls=x2mx3*x2mx3+y2my3*y2my3+z2mz3*z2mz3 
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r2s=xlmx3*xlmx3+ylmy3*ylmy3+zlmz3*zlmz3 

r3s=x2mxl*x2mxl+y2myl*y2myl+z2mzl*z2inzl 

c  compute  the  height  to  base  ratios. 
area2=area+area 
htbl=area2/rls 
htb2=area2/r2s 
htb3=area2/r3s 
htbmin=aminl (htbl ,htb2 ,htb3) 
tarea=tarea+area 
if ( ar ea . gt .  are  amx  )  then 
mxarea=if ace 
areamx=area 
endif 

if  (area .  It .  areamn)  then 
mnarea=if ace 
areamn=area 
endif 

if  (ht  bmin .  It .  rat  io  )  then 
mnrtio=iface 
ratio=htbmin 
endif 

40  continue 

avarea=tarea/nfaces 

write(3,110) 

110  f ormat(/25x,  'body  parameter  list'/) 
write (3 ,111)  nnodes ,nedges ,nf aces ,nunknb 

111  f ormat(l Ox, * number  of  vertices^ , i4,/10x, ’number  of  edges=' ,i4,/10 
>x, 'number  of  faces=' ,i4,/10x, ' number  of  edges  including  multiplici 
>ty=' , i4) 

write(3,205) 

205  format(/25x,  '  modeling  parameter  list  (meters)'/) 
write(3,206)  tarea 

206  format ( lOx , '  surface  area  of  the  scatterer=> ,el2. 5, lx, 'sq. meters' ) 
write (3 , 209 )  avedge , mxedge , edgemx , mnedge , edgemn 

209  format (lOx, 'average  edge  length^' , lel2 . 5 , lx, 'meters1 , 

$/10x, 'maximum  edge  length(edge  no. '  ,i3, 1 )=' ,el2.5,lx, 'meters' , 
$/10x, 'minimum  edge  length(edge  no. ' ,i3, ' )=' ,el2.5,lx, 'meters') 
write (3, 210)  avarea,mxarea,areamx,mnarea, areamn 

210  format (lOx, 'average  face  area  =' ,el2. 5, lx, 'sq. meters' ,/10x, 

$ 'maximum  face  area  (face  no. ' , i4, lx, ')=' ,el2 .5, lx, 'sq. meters ' ,/ 
$10x, 'minimum  face  area  (face  no. ' ,i4,lx, ')=' ,el2.5,lx, 'sq. meters ' ) 
write(3,211)  mnrtio, ratio 

211  format(l Ox, 'minimum  face  height  to  base  ratio  (face  no.', 

$i4,lx, ' )=' ,ell.5) 

return 

end 

subroutine  f acedg(nf aces , nbound , if ace , iedgl , iedg2 , iedg3) 
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c  this  subroutine  returns  the  edges  of  face  number  if ace. 
integer  nbound(3,6000) 
iedgl=nbound(l , if ace) 
iedg2=nbound(2 , iface) 
i edg3  =nbound ( 3 , if ac  e ) 
return 
end 

subroutine  f acvtx(nconn,nedges ,  iel ,ie2 ,  ie3 ,nvl  ,nv2 ,nv3) 

c  this  subroutine  returns  the  vertices  of  a  given  face 
c  nvl  is  node  opposite  edge  iel. 

c  nv2  is  node  opposite  edge  ie2. 

c  nv3  is  node  opposite  edge  ie3. 

integer  nconn(3,6000) 

c  the  node  nvl  is  the  node  that  edges  2  and  3  have  in  common, 
c  the  node  nv3  is  the  other  node  on  edge  2. 

if (nconn(l,ie2) .eq.nconn(l,ie3) .or .nconn(l ,ie2) .eq.nconn(2,ie3)) 
$then 

nv  1  =nconn  ( 1 ,  i  e2  ) 
nv3=nconn(2 , ie2) 
else 

nvl=nconn(2 ,  ie2) 
nv3=nconn(l , ie2) 
endif 

c  the  node  nv2  is  the  node  that  edges  1  and  3  have  in  common. 

if  (nconn(l,iel)  ,eq.nconn(l,ie3)  ,or.nconn(l  ,iel) . eq.nconn(2,ie3) ) 
$then 

nv2=nconn(l ,iel) 
else 

nv2=nconn(2 ,  iel) 
endif 
return 
end 

subroutine  vtxcrd(datnod , nnodes , nl , n2 , n3 , xl , x2 , x3 , yl , y2 , y3 , zl , z2 , 
>z3) 

c  this  subroutine  gets  the  coordinates  of  the  vertices  of  a  face 
dimension  datnod(3,6000) 
xl=datnod(l,nl) 
yl=datnod(2,nl) 
zl=datnod(3,nl) 
x2=datnod(l,n2) 
y2=datnod(2,n2) 
z2=datnod(3,n2) 
x3=datnod(l,n3) 
y 3=dat nod ( 2 , n3 ) 
z3=datnod(3,n3) 
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return 

end 


Subroutines  GEOM,  FACETCK,  and  AXB  are  used  here.  Listings  appear  in  Ap¬ 
pendix  A. 
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APPENDIX  C:  PATCH-TO-ACAD  TRANSLATOR  CODE 


As  described  in  Section  5.2,  there  are  two  versions  of  the  PATCH-to-ACAD  trans¬ 
lator:  PTA  and  PTF.  PTA  translates  the  entire  patch  file  geometry  as  a  single  part. 
PTF  translates  each  facet  as  a  single  part.  Only  PTF  is  listed  here.  Therefore,  in  the 
case  of  PTF,  the  file  length  is  much  greater,  but  ACAD  is  able  to  manipulate  each 
facet.  Facet  checking  is  used  in  subroutine  GEOM. 


c  program  ptf.f 

c  (similar  to  pta.f  version  3  -  uses  INSCRIBED  CIRCLE  face  check) 
c 

c  "patch  to  facet"  translator 

c  ******  differs  from  pta.f  in  that  each  triangle  is  ********* 
c  ******  an  individual  part  ********* 

c 

c  this  program  reads  a  file  named  "in. patch"  then  reformats  the  data 
c  and  writes  it  to  the  file  named  "out .facet"  which  can  be  read  into 
c  acad  as  a  *. facet  file 

c  x,y,z  are  node  coordinates  (index  is  facet  number) 
c  datnod  and  nconn  are  same  as  in  patch. f 

dimension  node (3 ,6000) ,nbound(3,6000) ,np(6000) , indsum(6000) 
dimension  datnod (3, 6000) ,nconn(3, 6000) , istart (6000) 
character *80  title 
character*  19  duml 
character* 11  dum2 

c  verbose  mode:  iverb=0  displays  progress 
iverb=0 

c  *******************  read  patch  file  ******************************* 
open ( 2 , f il e= , in . pat  ch 9 , s t atus=  *  old 9 ) 
read(2,ll)  title 
11  format (a80) 

read(2,*)  nverts ,nedges 
do  151  nn=l, nverts 

151  read (2,*)  n, datnod (1 ,n) , datnod (2, n) , datnod (3, n) 
do  200  nn=l, nedges 

200  read(2,*)  n,nconn(l,n)  ,nconn(2,n) 

c  fill  array  nbound  with  edges  of  each  face 

call  geom (datnod, nconn, nedges , indsum, nbound , 

$  nfaces ,nunknb) 
do  201  nf=l, nfaces 
istart (nf )=nf 
istart (nf )=nf+l 

call  prntbnd(nconn, nbound, istart ,nf, nedges, nfaces , nfaces) 

201  continue 

c  get  body  parameters 

call  bodpar (datnod , nconn , nbound , nverts , nedges ,nf aces , nunknb ) 
c  find  the  verticies  of  each  face  and  put  in  array  node 
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do  600  n=l ,nf aces 

if (iverb.eq.O)  write(6,*)  'finding  vertices  of  face  '  ,n 
il=nbound(l,n) 
up ( 1 ) =nconn ( 1 , il ) 
np(2)=nconn(2,il) 
i2=nbound(2,n) 
np(3)=nconn(l, i2) 
up  (  4  )  =n  c  onn  (  2 ,  i  2  ) 
i3=nbound(3,n) 
np(5)=nconn(l,i3) 
np  (  6  )  =n  c  o  nn  (  2 ,  i  3  ) 
c  find  the  three  unique  points 
node(l,n)=np(l) 
node(2,n)=np(2) 
ieg=3 

do  603  ii=3,6 
npt=np(ii) 

if  ((npt.ne.np(l))  .and.  (npt  .ne.np(2)) )  then 
c  must  be  the  third  node 
node(3,n)=npt 
go  to  602 
endif 

603  continue 

602  continue 

if (iverb.eq.O)  write(6,*)  'vertices  are:  ’ , (node(j j ,n) , j j=l ,3) 
600  continue 

c  *******************  yrite  facet  file  ******************************* 
titled FACET  FILE  V3.0  SG4D  ' 

open(l , file= 'out .facet' ) 
c  ignoring  material  parameters 
write(l,i)  title 

c  number  of  parts  is  the  number  of  facets 
write(l,2)  nfaces 
do  20  n=l, nfaces 
write (1,13)  '  FFace  '  ,n 
13  format (a7,i4) 
write(l,4)  0 

c  each  facet  has  3  vertices 
write(l,4)  3 
do  40  m=l , 3 

write(l ,8)  datnod (1, node (m, n) )  ,datnod(2, node (m,n))  , 

&  datnod(3,node(m,n)) 

40  continue 

write (1,2)  1 
write(l,*)  'Tri  Sheet  0' 
write(l,7)  3, 1,0, 0,0, 0,0 
write(l ,9)  1,2,3 
20  continue 

1  format (a28) 
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2  format (i5) 

3  format (a80) 

4  format (il) 

5  format (i7) 

6  format(2(f 14. 6, lx) ,f 14.6) 

7  format (7i7) 

8  format (3 (f 8.2, lx)) 

9  format (2x,3i5) 

c  write  new  data  in  MATLAB  files 
open(12,file=,xpts  .mO 
open(13,file='ypts .m' ) 
open ( 14 , f ile=  * zpt s . m  > ) 
open(15,f ile=*  endl .m' ) 
open ( 16, f ile=,end2.m,) 

90  format  (i5) 

do  140  n=i,nverts 
write(12,100)  datnod(l,n) 
write (13, 100)  datnod(2,n) 
write (14, 100)  datnod(3,n) 

140  continue 

100  format (f 15.4) 

do  150  n=l,nedges 
write(15, 101)  nconn(l,n) 
write(16,101)  nconn(2,n) 

150  continue 

101  format(i5) 

998  continue 

stop 

end 

0========================:=====================:=:=:====:================ 

subroutine  prntbnd(nconn,nbound, istart , i,nedges ,nf aces ,nbodys) 
0==================:=======  =  ==============:===:=  ====:=======:============ 

c  this  subroutine  prints  the  edges  and  the  vertices  of  each  face, 
c  input: 

c  nconn  has  the  vertices  and  the  multiplicity  factor  for  each  edge, 
c  nbound  has  the  edges  for  each  face, 
c  istart  has  the  beginning  faces  for  each  body, 
c  i  is  the  present  body, 
c  nedges  is  the  total  number  of  edges, 

c  nfaces  is  the  total  number  of  faces, 

c  nbodys  is  the  total  number  of  bodys. 

integer  nconn (3, 7000) , nbound (3, 7000) , istart (7000) 
integer  nverts(7000,3) 
open(3,f ile=,facelist J ) 
open(9 ,f ile= }f acedat ' ) 
open(20,file=*nodel .mJ) 
open(21,f ile=,node2.mJ) 
open ( 22 , f ile=  * node3 . m  * ) 
c  loop  through  the  faces  of  this  body. 
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do  10  ilO=istart (i) ,istart (i+l)-l 

call  f acvtx (nconn ,nedges ,nbound( 1 ,  ilO) ,nbound(2,il0) , 
>nbound(3,ii0)  ,nvl,nv2,nv3) 

write(3,98)il0,nbound(i,il0),nbound(2,il0) ,nbound(3,il0) ,nvi ,nv2 
>,nv  3 

c  write  face  number  and  vertices  to  “facedat11 
write(9,*)il0,nvl,nv2,nv3 
write(20, 100)  nvl 
write (21, 100)  nv2 
write (22, 100)  nv3 
nverts(i!0, l)=nvl 
nverts(il0,2)=nv2 
nv  ert  s ( il 0 , 3 ) =n v3 
10  continue 

98  format  (lac,  'face'  ,i5,  *  has  edges1, 3i5,'  with  vertices1 ,3i5) 

100  format (i6) 
return 
end 

subroutine  bodpar (datnod, nconn, nbound, nnodes ,nedges,nf aces, nunknb) 

c  this  subroutine  computes  the  body  parameters 
c  input: 

c  datnod(i, j)  i=l,2,3  are  the  x,y,z  coordinates  of  the  jth  node, 
c  nconn (3,  j )  :  edge  j  runs  from  node  nconn(l,j)  to  node  nconn(2,j) 
c  nconn (3, j)  is  the  multiplicity  of  the  jth  edge, 

c  nbound(i, j) :  contains  the  ith  edge  of  the  jth  face  i=l,2,3. 
c  nnodes  =  the  number  of  body  nodes . 

c  nedges  =  the  number  of  edges, 

c  nfaces  =  the  number  of  faces, 

c  nunknb  =  the  number  of  body  unknowns . 

c  output: 

c  avedge  =  the  average  edge  length(meters**2)  including  multiplicity, 
c  edgemx  =  the  maximum  edge  length (met ers ) . 

c  mxedge  =  the  edge  number  of  the  edge  with  length  edgemx. 

c  edgemn  =  the  minimum  edge  length  (met  ers)  . 

c  mnedge  =  the  edge  number  of  the  edge  with  length  edgemn. 

c  tarea  =  the  surface  area  of  the  scatter (meters**2) : for  thin 
c  structures  only  one  side  is  considered  in  the  surface  area, 

c  avarea  =  the  average  area  of  the  faces. 

c  mxarea  =  the  number  of  the  face  with  the  maximum  area(areamx) . 

c  mnarea  =  the  number  of  the  face  with  the  minimum  area(areamn)  . 

c  ratio  =  the  minimum  height  to  base  ratio  over  all  faces, 

c  mnrtio  =  the  face  number  that  has  a  height  to  base  ratio  of  'ratio1, 

dimension  datnod (3, 6000) 
integer  nconn(3 ,6000) ,nbound(3 ,6000) 

common/params/ avedge ,  edgemx,  mxedge ,  edg  emn,  mnedge ,  tarea,  avarea, 
$mxarea,  mnarea,  areamx,  areamn,  ratio ,  mnrtio 
c  ommon/ mchval/ valmax , valmin 
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c  save  /params/ 

c  the  following  line  is  a  statement  function. 

s ize (x  , y ,  z) =sqrt (x*x+y*y+z*z) 
c  initialization. 
sedgl=0 
edgemx=valmin 
edgemn=valmax 
do  20  ie=l,nedges 
mult=nconn(3 , ie) 
nl=nconn(l,ie) 
n2=nconn(2,ie) 
x=datnod ( 1 , n2 ) -datnod ( 1 , nl ) 
y =datnod ( 2 , n2 ) -datnod ( 2 , nl ) 
z=datnod ( 3 , n2 ) -datnod ( 3 , nl ) 
edgl=size (x , y , z) 
s  edgl=s  edgl+mult *  edgl 
if ( edgl . gt . edgemx) then 
edgemx=edgl 
mxedge^ie 
endif 

if  (  edgl .  It .  edgemn)  then 
edgemn=edgl 
mnedge=ie 
endif 

20  continue 

av  edg  e= s  edgl/nunknb 

c  compute  tarea,avarea,mnarea,areamn,mxarea,areamx, ratio,  and  mnrtio. 
ratio=valmax 
areamx=valmin 
ar  eamn= valmax 
tarea=0 . 

do  40  iface=l,nfaces 
isl=nbound(l , iface) 
is2=nbound(2 , iface) 
is3=nbound(3 , iface) 

call  f acvtx (nconn , nedges , is 1 , is2 , is3 , nvl , nv2 ,nv3 ) 

call  vtx crd( datnod ,nnodes , nvl ,nv2 , nv3 , xl , x2 ,x3 , y 1 , y2 , y3 , zl ,z2 , z3 

» 

xlmx3=xl-x3 

ylmy3=yl-y3 

zlmz3=zl-z3 

x2mx3=x2-x3 

y2my3=y2-y3 

z2mz3=z2-z3 

x2mxl=x2-xl 

y2myl=y2-yl 

z2mzl=z2-zl 

c  compute  area  of  face  by  taking  the  cross  product  of  two  edge  vectors. 
vx=y lmy 3*z2mz3-z lmz3*y 2my3 
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Vy=zlmz3*x2mx3~xlmx3*z2niz3 
vz=x lmx3 *y2my3-y lmy 3*x2mx3 
area= . 5* size ( vx , vy , vz) 

c  compute  the  square  of  the  lengths  of  each  side. 
rls=x2mx3*x2mx3+y2my3*y2my3+z2mz3*z2mz3 
r2s=xlmx3*xlmx3+ylmy3*ylmy3+zlmz3*zlmz3 
r3s=x2mxi*x2mxl+y2myl*y2myl+z2mzl*z2mzl 
c  compute  the  height  to  base  ratios. 
area2=area+area 
htbl=area2/rls 
htb2=area2/r2s 
htb3=area2/r3s 
htbrain=aminl(htbl ,htb2 ,htb3) 
tarea=tarea+area 
if ( ar ea . gt . ar eamx ) then 
mxarea=if ace 
areamx=area 
endif 

if  (area .  It .  areamn)  then 
mnarea=iface 
areamn=area 
endif 

if  (htbmin.lt  .ratio ) then 
mnrtio=iface 
ratio=htbmin 
endif 

40  continue 

avarea=tarea/nfaces 

write(3,110) 

110  format(/25x,  'body  parameter  list'/) 
write (3, 111)  nnodes  ,nedges  ,nfaces,nunknb 

111  format(10x, 'number  of  vertices=' ,i4,/10x, 'number  of  edges=' , i4,/10 
>x, 'number  of  f aces=' , i4, /10x, 'number  of  edges  including  multiplici 
>ty=' ,i4) 

write(3,205) 

205  format(/25x, 'modeling  parameter  list  (meters)'/) 
write (3 ,206)  tarea 

206  format (lOx, 'surf ace  area  of  the  scatterer=' ,el2.5,lx, 'sq. meters ') 
wr it e ( 3,209)  avedge , mxedge , edgemx , mnedge , edgemn 

209  format(10x, 'average  edge  length=' , lel2. 5, lx, 'meters' , 

$/10x, 'maximum  edge  length(edge  no. ' ,i3, ')=' ,el2.5,lx, 'meters' , 

$/i Ox, 'minimum  edge  length(edge  no. ' ,i3, ')=' ,el2.5,lx, 'meters') 
write (3, 210)  avar ea, mxarea,ar eamx, mnarea, areamn 

210  formatClOx, 'average  face  area  =' ,el2. 5, lx, 'sq. meters' ,/10x, 

$ 'maximum  face  area  (face  no. ' , i4, lx, ' )=' ,e!2 . 5, lx, 'sq. meters ' , / 
$10x, 'minimum  face  area  (face  no. ' ,i4,lx, ')=' ,el2.5,lx, 'sq. meters') 
write(3,211)  mnrtio, ratio 

211  format(10x, 'minimum  face  height  to  base  ratio  (face  no.', 

$i4,lx, ')=' ,ell .5) 
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return 

end 

subroutine  facedg(nfaces ,nbound,if ace , iedgl , iedg2, iedg3) 

c  this  subroutine  returns  the  edges  of  face  number  if ace. 
integer  nbound(3t7000) 
iedgl=nbound( 1 , if ace) 
iedg2=nbound(2 , if ace) 
iedg3=nbound(3 , if ace) 
return 
end 

subroutine  f acvtx(nconn,nedges , iel , ie2 , ie3 , nvl ,nv2 ,nv3) 

c  this  subroutine  returns  the  vertices  of  a  given  face 
c  nvl  is  node  opposite  edge  iel. 

c  nv2  is  node  opposite  edge  ie2. 

c  nv3  is  node  opposite  edge  ie3. 

integer  nconn(3,7000) 

c  the  node  nvl  is  the  node  that  edges  2  and  3  have  in  common, 
c  the  node  nv3  is  the  other  node  on  edge  2. 

if  (nconn(l,ie2)  .  eq.nconn(l,ie3)  . or  ,nconn(l , ie2)  .  eq.nconn(2, ie3) ) 
$then 

nvl=ncorm(l  ,ie2) 
nv3=nconn(2 , ie2) 
else 

nvl=nconn(2 ,  ie2) 
nv3=nconn ( 1 , ie2) 
endif 

c  the  node  nv2  is  the  node  that  edges  1  and  3  have  in  common. 

if  (nconn(l,  iel)  .  eq.nconn(l ,ie3)  .or  .nconn(l,iel)  .eq.nconn(2,ie3)) 
$then 

nv2 =nconn ( 1 , i e 1 ) 
else 

nv2=nconn ( 2 , i e 1 ) 
endif 
return 
end 

subroutine  vtxcrd(datnod,nnodes,nl,n2,n3,xl,x2,x3,yl,y2,y3,zl,z2, 
>z3) 

c  this  subroutine  gets  the  coordinates  of  the  vertices  of  a  face 
dimension  datnod(3, 6000) 
xl=datnod(l,nl) 
yl=datnod(2,nl) 
zl=datnod(3 ,nl) 
x2=datnod(l,n2) 
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y2=datnod(2,n2) 

z2=datnod(3,n2) 

x3=datnod(l,n3) 

y3=datnod(2,n3) 

z3=datnod(3,n3) 

return 

end 


Subroutines  GEOM,  FACETCK,  and  AXB  are  used  here.  Listings  appear  in  Ap¬ 
pendix  A. 
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APPENDIX  D:  PATCH-TO-NEC  TRANSLATOR  CODE 


The  PATCH-toNEC  translator  allows  the  user  to  convert  a  PATCH  input  file 
into  a  format  that  can  be  read  by  NEC.  CAUTION:  The  correspondence  is  only  ap¬ 
proximate  and  some  errors  may  result  as  noted  in  Section  5.3 

c  program  ptn.f 
c 

c  "patch  to  NEC"  translator 
c 

c  this  program  reads  a  file  named  "in. patch"  then  reformats  the  data 
c  and  writes  it  to  the  file  named  "out. nee"  which  can  be  read  by  NEC. 
c  x,y,z  axe  node  coordinates  (index  is  facet  number) 
c  datnod  and  nconn  are  same  as  in  patch. f 

c  if  iflag=0  edges  in  the  xy  plane  (i.e.,  z  <  eps)  are  deleted 
c  (omitted  in  the  NEC  GW  list) 

dimension  datnod (3, 8000) ,nconn(3,8000) 
character*80  title 
iflag=0 
eps=l . e-4 

c  *******************  read  patch  file  ******************************* 
open^jfile^in.patch*  ,status=Jold> ) 
read(2,12)  title 
12  format (a80) 

r ead ( 2 , *)  nvert s , nedges 
do  151  nn=l,nverts 

151  read(2,*)  n,datnod(l ,n) , datnod (2, n) ,datnod(3,n) 

do  200  nn=l, nedges 

200  read(2,*)  n,nconn(l,n)  ,nconn(2,n) 

c  *******************  write  NEC  file  ******************************* 
c  patch  edge  index  becomes  tag  number 
c  radius  of  each  wire  is  set  to  1/10  of  its  length 
open (l,file=J out .nee*) 
do  10  ne=l, nedges 
c  coordinates  of  first  end 
nl=nconn(l,ne) 
xl=datnod(l,nl) 
yl=datnod(2,nl) 
zl=datnod(3,nl) 
c  coordinates  of  second  end 
n2=nconn(2,ne) 
x2=datnod(l,n2) 
y2=datnod ( 2 , n2 ) 
z2=datnod(3,n2) 

c  if  iflag=0  and  this  segment  lies  in  the  xy  plane  omit  it  from 
c  the  GW  list.  Note  that  the  edge  list  is  not  compacted,  only 
c  they  axe  omitted  in  the  write 
if (iflag.eq.0)  then 
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if ((zl.lt.eps) .and. (z2.lt.eps))  go  to  20 
endif 

segl=sqrt((x2“Xl)**2+(y2-yl)**2+(z2--zl)**2) 

c  set  radius  for  this  segment 
r0=segl/10. 
if (rO.gt .0.2)  r0=0 . 2 
if (r0.lt .0.02)  r0=0.02 

write (1,300)  'GW' ,ne,l,xl,yl,zl,x2,y2,z2,r0 
20  continue 
10  continue 

300  format (a2, lx,2(i4, lx) ,6(f7.2,2x) ,f5.2) 
write(l ,400)  'GEM, 2 
write (1,400)  'GN',1 
write(l ,400)  JE1P 
400  format (a2,2x,il,2x,il) 


APPENDIX  E:  PATCH  INPUT  CHECKING  CODE  (KNIT) 

KNIT  checks  a  PATCH  input  file  for  duplicate  edges  and  nodes. 


c  program  knit.f  (final  version  4:  12/10/95) 
c 

c  removes  duplicate  nodes  in  a  inpatch  file  &  NONEDGES 
c  finds  lines  that  cross  or  are  parallel 
c  (an  edge  that  has  the  same  node  at  both  ends,  this 
c  allows  triagles  to  be  built  from  quadralaterals . ) 
c  edge  and  node  tags  do  not  have  to  be  ordered  in  , inknit } 
c  edge  numbers  out  will  be  different  than  edge  numbers  in 
c  datnod  is  the  same  as  in  patch.f 

dimension  newl(5000) ,new2(5000) , datnod (3, 5000) ,nnord(5000) 
dimension  tmpdat (3,5000) ,icount(5000) ,neord(5000) 
dimension  ivmin(5000) , ivtx(5000 , 2000) , iskip(5000) , indx(5000) 
dimension  nedge(SOOO) ,node(5000) ,nflag(5000) ,mf lag(5000) 
integer  tmpnod(5000) , end2 (5000) , endl (5000) 
character*80  title 
character* 8  fin,  fout 

data  nf lag/5000*0/ , mf lag/5000*0/ , iscl/1/ , iverb/0/ 
c  distances  less  than  eps  are  considered  the  same 
eps=l . e-3 

c  iverb=0:  verbose  mode  —  progress  displayed 
c  iscl=0:  rescale  data 
iverb=l 
xmin=l . e6 
ymin= 1 . e6 
zmin=l . e6 
xmax=-l .  e6 
ymax=-l.e6 
zmax=-l  ,e6 

c  read  the  inpatch  file  named  inknit 

***3^************************************************************** 

write(6,*) 

&  } enter  file  name  to  remove  duplicate  edges  and  nodes :* 
read(5,l)  fin 

1  format (a8) 
write(6,*) 

&  '  enter  output  file  name  (<9  char  &  different  from  input  name) ' 
read(5,l)  fout 

open( 1 , f ile=f in , status=  '  old* ) 
open(2,f ile=fout) 
read(l,2)  title 

2  format (a80) 
read(l,*)  nverts ,nedges 

write (6,*)  'nverts ,nedges=' , nverts ,nedges 
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nodes=nverts 
do  10  n=l,nverts 
read(l ,*)  node(n) ,xx,yy ,zz 
xmax=amax 1 ( xmax , xx ) 
ymax=amaxl  (yin  ax ,  yy) 
zmax=amax 1 (zmax , zz ) 
xmin=aminl (xmin , xx ) 
ymin=aminl  (ymin ,  yy  ) 
zmin= amin 1 ( zmin , z z ) 
tmpdat (l,n)=xx 
tmpdat (2, n)=yy 
tmpdat (3,n)=zz 
datnod ( 1 , n) =tmpdat ( 1 , n) 
datnod(2,n)=tmpdat(2,n) 
datnod (3 , n) - tmpdat (3  ,n) 
tmpnod (n) =node (n) 

c  if  the  node  has  been  read  change  the  flag  from  0  to  1 
if (nf lagCnode (n) ) . eq. 0)  nf lag (node (n))=l 

10  continue 

do  20  n=l,nedges 

read(l , *)  nedge(n) ,endl(n) ,end2(n) 

newl(n)=endl(n) 

new2(n)=end2(n) 

c  if  the  edge  has  been  read  change  the  flag  from  0  to  1 
if (mflag(nedge(n)) .eq.O)  mflag(nedge(n))=l 

20  continue 

c  check  to  see  if  any  nodes  have  not  been  read 
do  30  n=l,nverts 
if (nflag(n) .eq.O)  then 
write(6,*)  'node  index  skipped:  ',n 
endif 

30  continue 

c  check  to  see  if  any  edges  have  not  been  read 
do  35  n=l>nedges 
if (mflag(n) .eq.O)  then 
write(6,*)  'edge  index  skipped:  ' ,n 
endif 

35  continue 

write (6,*)  'finished  reading  file  inknit' 

c  order  node  numbers  from  lowest  to  highest  in  case  they  have  not 

c  been  defined  this  way 
do  573  nl=l, nodes 
do  570  n2=l, nodes 

c  find  edge  nl  and  save 

if (node(n2) .eq.nl)  then 
nnord(nl)=nl 

datnod ( 1 , nl ) =tmpdat ( 1 , n2 ) 
datnod(2,nl)=tmpdat(2,n2) 
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datnod  (  3 ,  nl )  =tmpdat  (3 ,  n2  ) 
go  to  572 
endif 

570  continue 

572  continue 

573  continue 

do  575  n=l, nodes 

tmpdat ( 1 , n) =datnod ( 1 , n ) 
tmpdat ( 2 ,  n) =datnod ( 2 , n) 
tmpdat (  3 ,  n) =datnod ( 3 , n) 
node (n) =nnord (n) 

575  continue 

write(6,*)  'nodes  re-ordered  from  low  to  high* 
c  find  duplicate  vertices  and  save  their  indices 
c  icount(i)  counts  the  number  of  times  vertex  i  occurs 
do  75  ivl=i,nverts 
x 1 = tmpdat ( 1 , i v 1 ) 
y l=tmpdat ( 2 , iv 1 ) 
z 1= tmpdat (  3 ,  i  v  1 ) 
icount(ivl)=l 

c  keep  track  of  vertex  number  ivl  and  its  duplicates 
c  ivtx(node  number,  occurance  number,  duplicate  index) 
ivtx(ivi , icount (ivl) )=ivl 
do  73  iv2=l,nverts 
if (ivl.ne. iv2)  then 
x2= tmpdat ( 1 , iv2) 
y2=tmpdat (2 , iv2) 
z2=tmpdat (3 , iv2) 
dx=abs(xi-x2) 
dy=abs(yl-y2) 
dz=abs(zl-z2) 

if  ((dx.lt  .eps)  .and.  (dy.lt.eps)  .and.  (dz.lt .  eps))  then 
icount (ivl) =icount (ivl) +1 
ivtx(ivl ,  icount  (  ivl  ))=iv2 
endif 
endif 

73  continue 

75  continue 

if (iverb.eq.O)  write (6,*)  'array  ivtx  filled' 
c  for  each  vertex  find  the  smallest  index 
do  78  iv=l,nverts 
ivmin ( i v ) =nvert  s + 1 
do  78  ii=l, icount (iv) 

ivmin ( iv) =min( ivmin ( iv ) , i vt x ( iv , ii ) ) 

78  continue 

if (iverb. eq.O)  write(6,*)  'found  smallest  index' 
c  find  duplicate  vertices 
irem=0 

do  76  ii=l,nverts 
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c  if  minimum  index  is  .It.  ii  this  is  a  duplicate 
if (ivmin(ii) .It .ii)  then 
irem=irem+l 
iskip(irem)=ii 
endif 

76  continue 

if (iverb. eq. 0)  write(6,*)  'duplicate  vertices  tagged' 
c  order  the  indices  to  be  removed  from  lowest  to  highest 

if (iverb. eq.O)  write(6,*)  'start  to  re-order  indices' 
idx=nverts 
do  85  il=l,irem 
do  83  i2=l,irem-il+l 
indx ( i 1 ) =min ( idx , iskip ( i2 ) ) 

83  continue 

idx=indx(il)+l 
85  continue 

if (iverb. eq.O)  then 
write(6,*)  'skip','  ', 'ordered' 
do  190  ix=l,irem 

190  write(6,*)  iskip(ix),'  ',indx(ix) 

endif 

c  fill  array  tmpnod  with  original  node  values 
do  196  nn=l,nverts 
tmpnod (nn) =node (nn) 

196  continue 

c  set  all  node  indices  to  minimum  value 
do  95  nn=l,n verts 
iv=tmpnod(nn) 
node (nn)=ivmin ( iv) 

95  continue 

c  refill  tmpnod  with  minimum  node  values 
do  96  nn=l,nverts 

c  set  all  endpoints  to  minimum  node  values 
do  91  kk=l,nedges 

if (endl(kk) .eq. tmpnod (nn))  newl (kk) =node (nn) 
if ( end2 (kk) . eq . tmpnod (nn) )  new2 (kk) =node (nn) 

91  continue 

96  continue 

do  192  kk=l,nedges 
endl (kk) =newl (kk) 
end2 ( kk ) =ne w2 (kk ) 

192  continue 

c  remove  duplicate  nodes  and  shift  remaining  nodes  down  one  for 
c  each  previous  node  removed 
do  79  ir=l,irem 
do  79  iv=indx(ir) ,nverts-ir 
datnod ( 1 , i v ) =datnod ( 1 , iv+ 1 ) 
dat nod ( 2 , i v ) =dat nod ( 2 , iv+ 1 ) 
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dat nod ( 3 , i v ) =dat nod ( 3 , i v+ 1 ) 

79  continue 

c  check  each  node  to  see  how  many  previous  nodes  have  been  removed 
c  (to  determine  the  number  of  steps  to  decrement  the  index) 
do  90  nn=l, averts 
idx=tmpnod (nn ) 

c  count  vertices  .It.  idx  that  have  been  removed 

c  nodes  greater  than  iskip(ir)  drop  one;  those  less  than  iskip(ir)  stay 
iter=0 

do  93  ir=l,irem 
if (idx.eq. indx(ir))  go  to  90 
do  93  ii=l,idx 

if ( ii . eq . indx ( ir) )  it er=iter+ 1 
93  continue 

it=it+i 

node ( i t ) = tmpnod ( nn ) - i t er 

c  find  all  edges  with  this  node  and  also  shift  them  down 
c  note  endl  &  end2  have  original  node  numbers;  newl  &  new2  have  new 
c  node  numbers 

do  179  kk=l,nedges 

if (endl (kk) . eq. tmpnod (nn))  newl(kk)=node(it) 
if (end2(kk) . eq. tmpnod (nn))  new2(kk)=node(it) 

179  continue 

90  continue 

nv er  t  s  =n  v  er  t  s  -  ir  em 

write(6,*)  'final  number  of  vertices  =', averts 
do  28  i=l,nedges 
endl(i)=newl(i) 

28  end2 ( i ) =new2 ( i) 

******************************************************************* 
c  rescale  data  if  desired 

******************************************************************* 
if (iscl.eq.O)  then 

write(6,*)  'rescale  data?  (0=yes/l=no) ' 
read(5,+)  ans 
if(ans.eq.O)  then 
write(6,*)  ' xmax , xmin= ' , xmax , xmin 
write (6,*)  'ymax,ymin=' ,ymax,ymin 
write(6,*)  'zmax,zmin=' ,zmax,zmin 
write (6,*)  'enter  scale  factor' 
read(5,*)  fac 
do  50  n=l, averts 
datnod(l,n)=fac*datnod(l ,n) 
datnod(2,n)=fac*datnod(2,n) 
dat nod ( 3 , n) =f  ac*dat nod ( 3 , n ) 

50  continue 
endif 
endif 

******************************************************************* 
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c  remove  duplicate  edges  and  use  new  node  numbers 
******************************************************************* 
c  order  edge  numbers  from  lowest  to  highest  in  case  they  have  not 
c  been  defined  this  way 
do  473  nl=l,nedges 
do  470  n2=l,nedges 
c  find  edge  nl  and  save 

if (nedge(n2). eq.nl)  then 
neord(nl)=nl 
v  newl (nl)=endl (n2) 
ne  w2  (n  1 ) = end2  (n2  ) 
go  to  472 
endif 

470  continue 

472  continue 

473  continue 

do  471  nn=l,nedges 
if (iverb.eq.O)  then 

write(6,*)  ' nn, nedge ,neord= *  ,nn,nedge(nn) ,neord(nn) 
endif 

nedge (nn) =neord (nn) 
endl (nn) =newl (nn) 
end2 (nn) =new2 (nn) 

471  continue 

write (6,*)  *  edges  re-ordered  from  low  to  high' 
********************* ********************************************** 
inon=0 

do  300  n=l,nedges 
if ( endl (n) . eq . end2 (n) )  then 
inon=inon+l 

c  nonedge  found  remove  edge  index  and  slide  all  edges  down  one 
do  310  nn=n,nedges+l-inon 
endl(nn)=endl(nn+l) 
end2  ( nn  ) = end2  (nn+ 1 ) 
newl(nn)=newl(nn+l) 
new2(nn)=new2(nn+l) 
nedge (nn) =nn 
neord(nn)=nn 
310  continue 

endif 

300  continue 

nedge s =n edge s-inon 

write(6,*)  'number  of  nonedges  found  was  ',inon 
if (iverb. eq.0)  then 

write (6,*)  ' start  to  generate  edge  connection  list* 
endif 

c  search  for  duplicate  edges 

c  for  each  edge  (n=l,nedges)  check  the  vertices  to  see  if  they  have 
c  been  previously  assigned  an  edge  number,  save  the  unique  set  of 
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c  endpoints  in  endl  and  end2.  length  will  be  ncount. 
do  175  iel=l,nedges 
icount (iel)=l 
npl=endl(iel) 
np2=end2(iel) 

c  keep  track  of  vertex  number  ivl  and  its  duplicates 
c  ivtx(node  number,  occurance  number,  duplicate  index) 
ivtx(iel, icount (iel))=iel 
do  173  ie2=l,nedges 
if (iel.ne.ie2)  then 
mpl=endl(ie2) 
mp2=end2(ie2) 

if (((npl.eq.mpl) .and. (np2. eq.mp2) ) .or. 

&  (Cnp2.eq.mpl) .and. (npl . eq.mp2) ) )  then 

c  .or. (npl .eq.np2))  then 

icount (iel)=icount(iel)+l 
ivtx(iel , icount (iel))=ie2 

endif 

endif 

173  continue 

175  continue 

c  for  each  edge  find  the  smallest  index 
do  178  ie=l,nedges 
ivmin(ie)=nedges+l 
do  178  ii=l, icount (ie) 

ivmin(ie)=min(ivmin(ie) , ivtx(ie, ii)) 

178  continue 
c  find  duplicate  edges 
irem=0 

do  176  ii=l,nedges 

c  if  minimum  index  is  .It.  ii  this  is  a  duplicate 
if (ivmin(ii) .It .ii)  then 
irem=irem+l 
iskip(irem)=ii 
endif 

176  continue 

c  order  the  edges  to  be  removed  from  lowest  to  highest 
write(6,*)  'number  of  duplicate  edges=',irem 
if (iverb.eq.O)  write (6,*)  're-ordering  edges' 
idx=n edges 
do  185  il=l,irem 
do  183  i2=l,irem-il+l 
indx ( i 1 ) =min ( idx , iskip ( i2 ) ) 

183  continue 

idx=indx(il)+l 
185  continue 

if (iverb. eq.O)  then 
write(6,*)  'skip','  'ordered' 
do  191  ix=l,irem 
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191 


' , indx(ix) 


write(6,*)  iskip(ix) , ' 
endif 

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
c  remove  duplicate  edges  and  shift  remaining  edges  down  one  for 
c  each  previous  edge  removed 
do  279  ir=l,irem 
do  279  ie=indx(ir) ,nedges-ir 
newl(ie)=newl(ie+l) 
new2(ie)=new2(ie+l) 
c  neord(ie)=neord(ie+l) 

279  continue 

nedges=nedges-irem 

write(6,*)  'final  number  of  edges  =',nedges 
if (iverb. eq.O)  write(6,*)  'edge  connection  list  generated' 

******  THIS  PORTION  OF  THE  CODE  HAS  NOT  BEEN  VALIDATED  ********* 
c  find  lines  that  cross  or  overlap 
c  write(6,*) 

C  &  >  ************************************************************  > 
c  do  909  na=l,nedges 

c  ne=nedge(na) 

c  nl=endl(ne) 

c  n2=end2(ne) 

c  xl=datnod(l,nl) 

c  yl=datnod(2,nl) 

c  zl=datnod(3,nl) 

c  x2=datnod(l,n2) 

c  y 2=datnod ( 2 , n2 ) 

c  z2=datnod(3,n2) 

c  do  908  mb=l,nedges 

c  me=nedge(mb) 

c  dont  check  if  same  edge 
c  if(me.le.ne)  go  to  907 

c  ml=endl(me) 

c  m2=end2(me) 

c  at  this  point  it  has  been  assumed  that  all  duplicate  edges  have 
c  been  removed,  don't  want  to  consider  lines  that  have  a  common  end  point 
c  if ((nl.eq.ml) .or. (n2.eq.m2) .or. 

c  &  (nl . eq.m2) .or . (n2. eq.ml))  go  to  907 
c  x3=datnod(l,ml) 

c  y3=datnod(2,ml) 

c  z3=datnod(3,ml) 

c  x4=datnod(l ,m2) 

c  y4=datnod ( 2 , m2 ) 

c  z4=datnod (3 , m2 ) 

c  write(6,*)  'calling  intersect  for  edges  ',me,ne 

c  call  intersect(xl,yl,zl,x2,y2,z2,x3,y3,z3, 

c  &  x4 , y4 , z4 , if lag , ipar ) 

c  if (ipar .eq.O) 

c  &  write(6,*)  'WARNING  -  edges  overlap:  ',me,ne 
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c  if (iflag.eq.O) 

c  ft  write(6,*)  'WARNING  -  edges  cross:  ',me,ne 

907  continue 

908  continue 

909  continue 
write(6,*) 

&  >  ^**************  +  *********************J|e*******J|este************  > 

**************************************************  j:**************** 
c  write  reformated  data  to  file 

open(2,f ile=fout) 
write(2,2)  title 
write (2,*)  nverts,nedges 
do  151  n=l,nverts 

151  write (2,*)  node (n)  , datnod ( 1 , n )  , datnod ( 2 , n)  , datnod ( 3 , n) 
do  200  n=l,nedges 
nedge (n) =neord (n) 

200  write (2,*)  nedge (n) ,newl(n) ,new2(n) 
c  write  new  data  in  MATLAB  files 
open(12,f ile='xpts .m' ) 
open ( 13 , f ile=  * ypt s . m  ' ) 
open(14,f ile=,zpts .m') 
open(15,f ile='endl .m') 
open(16,file='end2.m' ) 

98  format (i5) 

do  140  n=l,nverts 
write(12, 100)  datnod(l,n) 
write(13, 100)  datnod(2,n) 
write(14, 100)  datnod(3,n) 

140  continue 

100  f ormat (f 15 .4) 

do  150  n=l,nedges 
write(15,101)  newl(n) 
write(16,101)  new2(n) 

150  continue 

101  format (i5) 

998  continue 

stop 

end 

subrout ine  int er se  ct (xl , y 1 , z 1 , x2 , y 2 , z2 , x3 , y3 , z3 , 

&  x4,y4,z4, if lag,ipar) 

c  subroutine  to  find  the  intersection  of  two  lines 
c  endpoints  of  line  1:  (xl,yl,zl)  and  (x2,y2,z2) 
c  endpoints  of  line  2:  (x3,y3,z3)  and  (x4,y4,z4) 
c  if  there  is  an  intersection  iflag=0 
c  ipar=0  denotes  overlapping  lines 
tol=l . e-5 
if lag=l 
ipar=l 
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al=x2-xl 

bi=y2-yl 

ci=z2-zl 

a2=x4-x3 

b2=y4-y3 

c2=z4-z3 

vl=bl*c2-b2*cl 

v2=cl*a2-c2*al 

v3=a2*bl-al*b2 

dl=x3-xl 

d2=y3-yi 

d3=z3-zl 

rdotv2=dl*a2+d2*b2+d3*c2 

rdotvl=dl*al+d2*bl+d3*cl 

vldotv2=al*a2+bl*b2+cl*c2 

v2dotv2=a2*a2+b2*b2+c2*c2 

vlxv2sq=vl**2+v2**2+v3**2 

if (vlxv2sq.lt .tol)  vlxv2sq=0. 

c  check  to  see  if  the  lines  are  parallel;  set  ipar=0  if  cross 
c  product  is  zero 

if (vlxv2sq. eq. 0 . )  then 

c  if  these  lines  are  parallel  see  if  they  overlap,  note  it  is 
c  assumed  that  this  subroutine  is  not  called  for  the  same  edge 
c  *******  case  1:  none  of  the  coeffs  are  zero  ******* 
if ((abs(al) .gt .tol) .and. (abs(bl) .gt .tol) . and. 

&  (abs(cl) .gt .tol))  then 
c  write(6,*)  'easel* 

p=(x3-xl)/al 
q=(y3-yl)/bl 
r=(z3-zl)/cl 

if ( (abs(p-q) .gt .tol) .and. (abs(q-r) .gt. tol) .and. 

&  (abs(p-r) .gt .tol))  ipar=0 

p=(x4-xl)/al 
q=(y4-yl)/bl 
r=(z4-zl)/cl 

if ( (abs(p-q)  .gt .tol) .and. (abs(q-r) .gt .tol) .and. 

&  (abs (p-r) .gt .tol))  ipar=0 

endif 

c  *******  case  2:  one  of  the  coeffs  are  zero  ******* 
c  if  a 1=0  compare  y  and  z 

if ((abs(al) .It. tol) .and. (abs(bl) .gt .tol) .and. 

&  (abs (cl) .gt .tol))  then 
c  write(6,*)  *case2a* 

q=(y3-yl)/bl 
r=(z3-zl)/cl 

if (abs(r-q) .gt .tol)  ipar=0 

q=(y4-yl)/bl 

r=(z4-zl)/cl 

if (abs(r-q) .gt .tol)  ipar=0 
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endif 

c  if  bl-0  compare  x  and  z 

if ((abs(bl) .It .tol) .and. (abs(al) .gt.tol) .and. 

&  (abs (cl) .gt.tol))  then 
c  write(6,*)  ,case2b' 

p=(x3-xl)/al 
r=(z3-zl)/cl 

if (abs (p-r) .gt.tol)  ipar=0 

p=(x4-xl)/al 

r=(z4-zl)/cl 

if (abs (p-r) .gt.tol)  ipar=0 
endif 

c  if  cl=0  compare  x  and  y 

if ((abs(cl) .It .tol) .and. (abs(bl) .gt .tol) . and. 

&  (abs(al) .gt.tol))  then 
c  write(6,*)  ,case2c> 

p=(x3-xl)/al 
q=(y3-yl)/bl 

if (abs (p-q) .gt.tol)  ipar=0 

p=(x4-xl)/al 

q=(y4-yl)/bl 

if (abs (p-q) .gt.tol)  ipar=0 
endif 

c  *******  case  3:  two  of  the  coeffs  are  zero  ******* 
c  if  al=0  and  bl=0  line  is  parallel  to  z  axis 

if ((abs(bl) .It .tol) .and. (abs(al) .It .tol))  then 
if ((abs(xl-x3) .It .tol) .and. (abs(yl-y3) .It .tol))  ipar=0 
c  write(6,*)  'caseSa* 

endif 

c  if  al=0  and  cl=0  line  is  parallel  to  y  axis 

if ((abs(cl) .It .tol) .and. (abs(al) .It .tol))  then 
if ((abs(xl-x3) .It .tol) .and. (abs(zl-z3) .It .tol))  ipar=0 
c  write(6,*)  'caseSb' 

endif 

c  if  bl=0  and  cl=0  line  is  parallel  to  x  axis 

if  ((abs(cl)  .It  .tol)  .and.  (abs(bl)  .It. tol))  then 
if  ((abs(yi-y3)  .It  .tol)  .and.  (abs(zl-z3)  .It .tol))  ipar=0 
c  write (6,*)  'caseSc* 

endif 
endif 

100  continue 
c  find  intersections 

if (ipar .ne.O)  then 

tee=(rdotvl*v2dotv2-rdotv2*v2dotvl)/vlxv2sq 
if ((tee.gt .0. ) .and. (tee.lt . 1 . ))  then 
if lag=0 
xO=xl+tee*al 
yO=yl+tee*bl 
zO=zl+tee*cl 
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c  check: 


u0=x3+tee*a2 

v0=y3+tee*b2 

w0=z3+tee*c2 

if ((abs(xO-uO) .gt .tol) .or. (abs(yO-vO) ,gt . tol) . or. 
&  (abs(zO-wO) .gt .tol))  iflag=l 
endif 
endif 

200  return 
end 
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APPENDIX  F:  GEOMETRY  FILE  BUILDER  (BLDMAT) 

BLDMAT  generates  data  files  with  the  “.m”  extension  so  that  they  can  be  loaded 
into  MATLAB.  The  files  are  used  by  the  matlab  script  PLTPATCH  to  view  a  three- 
dimensional  plot  of  the  geometry.  BLDMAT  also  generates  and  edge  connection  list 
that  is  identical  to  that  of  PATCH. 


c  program  bldmatf ck . f  (DEVELOPMENT  version  2:  2/96) 

c  program  to  read  inpatch  data  (output  from  buildnS  or  mbuild) 
c  and  write  the  variables  to  *.m  files  for  use  by  MATLAB 
*********************************************************** 

parameter (mxunkn=5500, mxbdnd=2500 ,mxedgs=5500,mxface=3700 , 
$  mxfu=5500,mxdjbd=50,mxmult=3) 
*********************************************************** 
dimension  x(mxbdnd)  ,y(mxbdnd)  ,z(mxbdnd)  ,node(mxbdnd) 
dimension  nedge(mxedgs)  ,npl(mxedgs)  ,np2(mxedgs)  ,nbe(mxdjbd) 
dimension  datnod(3,mxbdnd)  ,nconn(3,mxedgs)  ,nbound(3,mxf ace) 
dimension  indsum(mxedgs)  ,istart(mxdjbd+l)  ,ipvt(mxfu) 
dimens ion  i edgf (mxmult + 1 , mxedgs ) 
character*50  title 
character*9  name in 
mult 1 =mxmult + 1 

write (6,*)  3 enter  input  file  name* 
read (5, 3)  namein 
3  format (a9) 

open(l  jfile^amein.statuss'oldO 
read(l,2)  title 
2  format (a50) 

read(l,*)  nvert,nedges 
nnodes=nvert 

write (6,*)  'nvert ,nedges=' ,nvert ,nedges 
do  10  n=l,nvert 

10  read(l,*)  node(n) ,x(n) ,y(n) ,z(n) 

do  20  n=l,nedges 

20  read(l,*)  nedge(n) ,npl(n) ,np2(n) 

c  write  new  data  in  MATLAB  files 
open(12,f  ile=*xpts  .mO 
open(13,f  ile=,ypts.mO 
open(14 , f ile= 3 zpt  s . m 3 ) 
open(15,f  ile=,endl.mO 
open(16 , f ile= 3 end2 .  m' ) 
open(3,file=,facef  ckO 
write (3,*)  namein 
90  format (i5) 

do  40  n=l,nvert 
write(12,100)  x(n) 
write(13,100)  y(n) 
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write(14,100)  z(n) 

40  continue 

100  format(fl5.4) 

do  50  n=i, nedges 
write(15, 101)  npl(n) 
write(16, 101)  np2(n) 

50  continue 

101  format (i5) 

c  fill  datnod  with  node  locations 
do  110  n=l,nvert 
datnod( 1 , n) =x (n) 
datnod( 2 , n) =y (n) 
datnod(3,n)=z(n) 

110  continue 

c  fill  nconn  with  edge  connections, 
do  120  n=l,nedges 
nconn(l ,n)=npl (n) 
nconn(2,n)=np2(n) 
nconn(3,n)=-l 
120  continue 

call  geom(datnod, nconn, nedges ,indsum, nbound, mxf ace, nf aces, 

$  mxbdnd,nunknb,mxedgs) 

call  curdir  (nconn ,  nbound ,  nf  aces  ,nedges ,  mxd  j  bd ,  ipvt ,  istart , 

$  nbodys,nbe) 

write(6,*)  ’number  of  disjoint  bodies  is  ’ ,nbodys 
do  210  i=l,nbodys 

call  prntbnd (nconn, nbound, istart , i ,nedges ,nf aces ,nbodys) 

210  continue 

call  edgf  ac  (nconn,  nedges,  nbound, nf  aces  ,iedgf, multi) 
c  call  edgep (nnodes , nedges , datnod , nconn , nunknb 
c  $  ,iedgf ,multl, nbound, nf aces) 
c  get  body  parameters 

call  bodpar (datnod, nconn, nbound, nnodes , nedges ,nf aces ,nbodys , 

$  nunknb) 
stop 
end 

subroutine  curdir (nconn, nbound, nf aces, nedges, mxd jbd,itree, is tart 
$nbodys  ,nbe) 

c  all  edges  in  the  first  disjoint  surface  are  numbered  consecutively 
c  starting  from  1.  the  edges  in  the  next  disjoint  surface  are  numbered 
c  consecutively,  starting  where  the  last  surface  left  off. 
c 

c  input: 

c  nconn(3, nedges) :  edge  j  runs  from  vertex  nconn(l,j)  to  vertex 
c  nconn(2,j).  nconn(3, j ) multiplicity  factor  of  the  edge, 
c  nbound ( 3, nf aces ): each  face  j  has  edges  nbound(i,j)  i=l,2,3 
c  j=l ,2, . . . ,nfaces . 


c  iif aces  equals  the  total  number  of  faces, 
c  nedges  equals  the  total  number  of  edges, 
c  mxdjbd  equals  the  maximum  number  of  expected  bodys. 
c 

c  output: 

c  istart (mxdjbd+1) : istart (i)=the  lowest  numbered  face  on  the  ith 
c  tree (disjoint  surface) 

c  istart (nbodys+l)=nfaces+l 

c  mxdjbd.ge.nbodys  or  routine  stops  and  prints  a  warning, 
c  itree(nfaces) 

c  itree(i)  i=l , . . . .istart (2)-l  =the  faces  on  the  first  tree, 

c  itree(i)  i=istart(j) , , . . , istart (j+l)-l,=the  faces  on  the  jth  tree, 

c  nbodys  equals  the  number  of  disjoint  surfaces. 

c  nbe(mxdjbd)  :nbe(i)  contains  the  number  of  boundry  edges  for  body  i. 
c 

integer  nconn (3, nedges) ,nbound(3,nf aces) ,itree(nfaces) , 

>  istart (raxdjbd+1) ,nbe(mxdjbd) 

c  set  present  tree  to  first  tree, 
c  ntree=total  number  of  faces  stored  in  the  tree, 
c  lnf=lowest  numbered  face  occuring  in  the  present  tree, 
c  conveniently  it  happens  that  itree(lnf )=lnf  so  lnf=lowest  index  i 
c  so  that  itree(i)  is  in  the  present  tree, 
do  40  i=l, mxdjbd 
nbe(i)=0 
40  continue 

nf acel=nf aces+1 

ntree=l 

Inf =1 

itree(lnf )=lnf 

istart (l)=lnf 

do  1  nbodys=l, mxdjbd 

c  add  the  number  of  boundry  edges  in  the  first  face  of  this  body  to 
c  nbe(nbody) . 

if (nconn(3 , nbound( 1 , Inf )). eq . 0) nbe (nbodys ) =nbe (nbodys ) +1 
if (nconn ( 3 , nbound ( 2 , Inf ) ) . eq . 0) nbe (nbodys ) =nbe (nbodys ) + 1 
if (nconn ( 3 , nbound ( 3 , Inf ) ) . eq . 0) nbe (nbodys ) =nbe (nbodys ) + 1 
c  search  for  a  face  that  may  be  added  to  the  present  tree. 

51  do  50  if ace=lnf+l,nfaces 
c  if  if ace  is  already  in  tree  continue  search, 
do  10  jtree=lnf .ntree 

if (iface.eq.itree(jtree))goto  50 
10  continue 

c  test  to  see  if  if ace  has  an  edge  in  common  with  present  tree. 
lowface=nfacel 

c  find  the  lowest  face  with  a  common  edge. 

do  20  jtree=lnf ,ntree 
c  test  for  a  common  edge. 

do  30  i=l ,3 
do  31  j=l, 3 


80 


if  (nbound  (  i ,  if  ace)  .  eq .  nbound  (  j ,  itr  e  e  (  j  tree  )  )  .  and . 
$  itree(jtree) .It . lowf ace)lowface=itree(j tree) 

31  continue 

30  continue 

20  continue 

if ( lowf  ac  e . ne . nf  ace 1 ) then 
c  common  edge  has  been  found. 


do  60  i=l ,3 
do  61  j=l,3 

if (nbound ( i , if  ace) . eq . nbound ( j , lowf  ace)) then 
c  incriment  number  of  faces  in  the  tree  and  add  if  ace  to  present  tree. 
ntree=ntree+l 
itree(ntree)=iface 

c  add  the  number  of  boundry  edges  in  this  face  to  nbe. 

if (nconn(3#nbound(l ,iface) ) . eq.O) 

>  nbe (nbody s ) =nbe (nbody s ) + 1 
if (nconn(3,nbound(2,iface)) .eq.O) 

>  nbe (nbody s)=nbe(nbodys)+l 
if ( nconn ( 3 , nbound (3, if ace)). eq.O) 

>  nbe (nbody s)=nbe(nbodys)+l 


goto  51 
endif 

61  continue 

60  continue 

endif 

50  continue 

lnf=ntree+l 
if (Inf .le.nf aces) then 
ist art (nbody s+l)=lnf 
c  initialize  new  tree. 

ntree=ntree+l 
itr ee (ntree) =lnf 
else 

goto  999 
endif 

1  continue 

write(3,99) 

99  f ormat( lx, * warning  in  curdir  mxdjbd  found  but  still  have  faces', 
$'  left 9 ) 


stop 

999  continue 

istart  (nbodys+l)=nf  aces+1 

return 

end 


subrout  ine  bodpar  (datnod ,  nconn ,  nbound ,  nnodes ,  nedges ,  nf  aces , 
$  nbody s,nunknb) 

c  input : 
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c  datnod(i,j)  i=l,2,3  are  the  x, y,z  coordinates  of  the  jth  node, 
c  nconn(3,j):  edge  j  runs  from  node  nconn(l,j)  to  node  nconn(2,j) 
c  nconn(3,j)  is  the  multiplicity  of  the  jth  edge, 

c  nbound(i,j):  contains  the  ith  edge  of  the  jth  face  i=l,2,3. 
c  nnodes  =  the  number  of  body  nodes, 

c  nedges  =  the  number  of  edges, 

c  nfaces  =  the  number  of  faces, 

c  nunknb  =  the  number  of  body  unknowns, 

c  output: 

c  avedge  =  the  average  edge  length(meters**2)  including  multiplicity, 
c  edgemx  =  the  maximum  edge  length (meters) . 
c  mxedge  =  the  edge  number  of  the  edge  with  length  edgemx. 
c  edgemn  =  the  minimum  edge  length  (meters ) . 
c  mnedge  =  the  edge  number  of  the  edge  with  length  edgemn. 
c  tarea  =  the  surface  area  of  the  scatter(meters**2) :f or  thin 
c  structures  only  one  side  is  considered  in  the  surface  area, 

c  avarea  =  the  average  area  of  the  faces. 

c  mxarea  =  the  number  of  the  face  with  the  maximum  area(areamx) . 

c  mnarea  =  the  number  of  the  face  with  the  minimum  area(areamn)  . 

c  ratio  =  the  minimum  height  to  base  ratio  over  all  faces, 

c  mnrtio  =  the  face  number  that  has  a  height  to  base  ratio  of  'ratio*, 

dimension  datnod( 3, nnodes) 
integer  nconn (3, nedges) ,nbound (3, nfaces) 
c  common/params/ avedge ,  edgemx , mxedge ,  edgemn,  mnedge ,  tarea ,  avarea, 

c  $mxar  ea ,  mnarea ,  areamx ,  areamn ,  rat  io ,  mnrt  io 

c  the  following  line  is  a  statement  function. 

size (x , y , z) =sqrt (x*x+y*y+z*z) 
c  initialization. 
sedgl=0 
valmax=l . e35 
valmin=-l . e35 
edgemx=valmin 
edgemn=valmax 
do  20  ie=l, nedges 
mult=nconn(3 , ie) 
nl=nconn(l,ie) 
n2=nconn(2,ie) 
x=datnod ( 1 , n2 ) -datnod ( 1 , nl ) 
y =datnod ( 2 , n2 ) -dat nod ( 2 , nl ) 
z=datnod ( 3 , n2 ) -datnod ( 3 , nl ) 
edgl=size(x,y,z) 
sedgl=s edgl +mult* edgl 
if ( edgl . gt . edgemx) then 
edgemx=edgl 
mxedge=ie 
endif 

if ( edgl . It . edgemn) then 
edgemn=edgl 
mnedge=ie 
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end  if 

20  continue 

av  edg  e = s  edgl/nunknb 

c  compute  t ax ea,avarea,mnarea,areamn,mxarea,areamx, ratio, and  mnrtio. 
ratio=valmax 
areamx=valmin 
ar eamn= valmax 
tarea=0 . 

do  40  iface=l ,nfaces 
isl=pbound( 1 , if ace) 
is2=nbound(2 , if ace) 
is3=nbound(3 , if ace) 

call  f acvtx (nconn , nedges , is 1 , is2 , is3 , nvl , nv2 , nv3 ) 

call  vtxcrd (datnod ,  nnodes ,  nvl ,  nv2 ,  nv3 ,  xl ,  x2 ,  x3 ,  y  1 ,  y2 ,  y3 ,  z  1 ,  z2 ,  z3 

» 

xlmx3=xl-x3 

ylmy3=yl-y3 

zlmz3=zl-z3 

x2mx3=x2-x3 

y2my3=y2-y3 

z2mz3=z2-z3 

x2mxl=x2-xl 

y2myl=y2~yl 

z2mzl=z2-zl 

c  compute  axea  of  face  by  taking  the  cross  product  of  two  edge  vectors. 
vx=y Imy3*z2mz3-z Imz3*y2my3 
vy=zlmz3*x2mx3-xlmx3*z2itiz3 
vz=xlmx3*y2my3-yimy3*x2mx3 
area=.5*size(vx, vy ,vz) 

c  compute  the  square  of  the  lengths  of  each  side. 
rls=x2mx3*x2mx3+y2my3*y2my3+z2mz3*z2mz3 
r2s=xlmx3*xlmx3+ylmy3*ylmy3+zlmz3*zlmz3 
r3 s =x2mx 1 *x2rax 1+y 2my 1 *y 2my 1 +z2mz 1 *z2mz 1 

c  compute  the  height  to  base  ratios. 
area2=area+area 
htbl=area2/rls 
htb2=area2/r2s 
htb3=area2/r3s 
htbmin=aminl (htbl , htb2 , htb3 ) 
tarea=tarea+area 
if ( area . gt . ar eamx ) then 
mxcirea=if  ace 
areamx=area 
endif 

if  (  area .  It .  ar  eamn)  then 
mnarea=iface 
areamn=area 
endif 

if  (htbmin .  It .  ratio)  then 
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mnrtio=iface 
ratio=htbmin 
endif 
continue 

avarea=tarea/nfaces 
write(3, 110) 

110  f ormat(/25x, 'body  parameter  list'/) 

write (3 ,111)  nnodes ,nedges ,nf  aces  ,nunknb,nbodys 

111  format(10x, 'number  of  vertices^' ,i5,/10x, 'number  of  edges= ' , i5, /10 
>x, 'number  of  faces=' , i5,/10x, 'number  of  edges  including  multiplici 
>ty=' ,i5,/10x, 'number  of  bodies=',i5) 

write (3, 205) 

205  format (/25x, 'modeling  parameter  list  (meters)'/) 
write(3,206)  tarea 

206  format(10x,  'surface  area  of  the  scattered'  ,el2. 5, lx,  'sq. meters') 
write (3 , 209)  avedge , mxedge , edgemx ,mnedge , edgemn 

209  format (lOx, 'average  edge  length=' ,le!2. 5, lx, 'meters' , 

{/10x, 'maximum  edge  length(edge  no . ' ,i5, ' )=' ,el2. 5, lx, 'meters' , 
$/10x, 'minimum  edge  length(edge  no. ' ,i5, ')=' ,el2.5,lx, 'meters') 
write (3, 210)  avarea,mxarea,areamx,mnarea,areamn 

210  format (10x, 'average  face  area  =' ,el2. 5, lx, 'sq. meters' ,/10x, 
{'maximum  face  area  (face  no. ' ,i5,lx, ')=' ,el2.5,lx, 'sq. meters' ,/ 
$10x, 'minimum  face  area  (face  no. ' ,i5,lx, ')  =  ' ,el2.5,lx, 'sq. meters') 

write(3,211)  mnrtio, ratio 

211  format(10x, 'minimum  face  height  to  base  ratio  (face  no.', 

{i5,lx, ')=' ,eli.5) 

if ((areamn.lt . l.e-10) .or. (ratio.lt . 1 .e-10) )  then 
write (6,*)  'TRIANGLES  WITH  ZERO  AREA' 
stop 
endif 
return 
end 


subroutine  f acvtx(nconn,nedges , iel , ie2 , ie3 ,nvl ,nv2 ,nv3) 

c  nvl  is  node  opposite  edge  iel. 

c  nv2  is  node  opposite  edge  ie2. 

c  nv3  is  node  opposite  edge  ie3. 

integer  nconn(3,nedges) 

c  the  node  nvl  is  the  node  that  edges  2  and  3  have  in  common, 
c  the  node  nv3  is  the  other  node  on  edge  2. 

if (nconn(l , ie2) . eq.nconn(l ,ie3) .or .nconn(l , ie2) . eq.nconn(2, ie3) ) 
{then 

nvl=nconn(l ,  ie2) 
nv3=nconn(2 ,  ie2) 
else 

nvl=nconn(2 ,  ie2) 
nv3=nconn(l  ,ie2) 
endif 
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c  the  node  nv2  is  the  node  that  edges  1  and  3  have  in  common. 

if  (nconn(l , iel)  . eq.nconn(l ,ie3)  .or .nconn(l ,iel)  .  eq.nconn(2,ie3)) 
$then 

nv2=nconn  ( 1 , ie 1 ) 
else 

nv2=nconn(2 ,  iel) 
endif 
return 
end 

subrout ine  vt x cr d ( dat nod , nnode  s , n 1 , n2 , n3 , x 1 , x2 , x3 , y 1 , y 2 , y 3 , z 1 , z2 , 
>z3) 

dimension  datnod(3,nnodes) 

xl=datnod(l,nl) 

yi=datnod(2,nl) 

zl=datnod(3 ,nl) 

x2=datnod(l,n2) 

y2=datnod(2,n2) 

z2=datnod(3,n2) 

x3=datnod(l,n3) 

y 3=datnod ( 2 , n3 ) 

z3=datnod ( 3 , n3 ) 

return 

end 

subrout ine  edgf ac (nconn , nedges , nbound , nf ace  s , iedgf , mult 1 ) 
c  input : 

c  edge  ie  runs  from  vertex  nconn(l,ie)  to  vertex  nconn(2,ie) 
c  and  has  multiplicity  nconn(3,ie). 
c  face  iface  has  edges  nbound( j ,if ace)  j=l,2,3 
c  multi  is  set  in  the  main  program  and  mult 1-1. ge. the 
c  maximum  multiplicity  of  any  edge, 
c  output : 

c  array  iedgf  for  an  edge  with  multiplicity  mult 
c  iedgf (l,ie)=the  lowest  numbered  face  connected  to  edge  ie. 
c  iedgf (2,ie)=the  next  lowest  numbered  face  connected  to  ie. 
c  . 

c  iegdf (mm,ie)=the  last  face  connected  to  ie. 
c  where  mm  is  the  number  of  faces  connected  to  edge  ie. 

integer  nconn(3, nedges) , nbound (3 ,nf aces ), iedgf (multi , nedges) 
c  initialize  the  array  iedgf. 
do  S  ie=l,nedges 
do  6  m=l, multi 
iedgf (m,ie)=0 
6  continue 

5  continue 
c  fill  array  iedgf. 
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do  100  ie=l, nedges 
multel=nconn(3 , ie)+l 
m=0 

do  50  if=l,nfaces 

if(m.ge.multel)go  to  100 

if (ie . eq. nbound( 1 , if ) . or . ie . eq. abound (2 , if) . or . ie . eq. 

$  nbound(3,if))then 
m=m+l 

iedgf (m,ie)=if 
endif 

50  continue 
100  continue 

write(3,*)  ’  7 
do  200  ie=l,nedges 
write(3,201)ie 

write(3,*)(iedgf (m,ie)  ,m=l,nconn(3,ie)+l) 

200  continue 

201  format (lx, 'edge*  ,i5,  *  is  attached  to  faces  1 ,) 
return 

end 

subrout  ine  edgep  (nnodes ,  nedge s ,  dat nod ,  nconn ,  nunknb 

>  ,iedgf , multi  ,nbound,nfaces) 

c  input: 

c  nnodes=the  number  of  body  nodes, 
c  nedges=the  number  of  edges. 

c  nunknb=the  number  of  body  unknowns  before  considering  the 
c  symmetry  plane  attachments. 

c  multl=the  maximum  allowed  multiplicity  for  an  edge  plus  1. 
c  nfaces=the  number  of  faces. 

c  datnod(i,n)=the  x,y,z  components (i=l ,2,3)  of  the  nth  node  n=l, nnodes. 
c  nconn(i,ie)  i=i,2,3:  edge  ie  (ie=l ,nedges)  runs  from  node  nconn(l,ie) 
c  to  nconn(2,ie)  and  has  multiplicity  nconn(3,ie) (before  any  symmetry 
c  plane  attachments  are  considered) . 

c  nbound(j ,  if  ace)  j=l,2,3  are  the  three  edges  attached  to  face 
c  number  if ace.  iface=l ,nf aces . 

c  iedgf(m,ie)  m=l, . . . ,nconn(3,ie)+l  contains  the  faces  attached  to  edge 
c  number  ie  (before  any  symmetry  plane  attachments  are  considered), 
c  ie=i, . . . ,nedges. 
c  output: 

c  for  each  edge  ie  that  is  connected  to  at  least  one  a  p.e.c.  symmetry  plane 
c  the  number  of  body  unknowns  (nunknb)  is  incremented  by  1. 
c  the  edge  vertex  connection  list  with  edge  multiplicities  is  outputted 
c  after  accounting  for  all  symmetry  plane  attachments. 

dimension  dat nod (3, nnodes) ,nconn(3,nedges) ,igndp(3) , 

>  iedgf (multi, nedges) ,nbound(3,nfaces) 
common/gplane/  ngndp ,  igndp 

logical  lfcpmc 
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save  /gplane/ 
if (ngndp . gt . 0 ) then 

c  if  the  maximum  distance  from  an  edge  to  a  symmetry  plane  is  less  than  or 

c  equal  to  edged,  the  edge  is  assumed  connected  to  that  symmetry  plane. 
edged=le-7 
do  100  ie=l,nedges 
nl=nconn(l,ie) 
n2=nconn(2,ie) 
xl=abs(datnod(l ,nl)) 
yl=abs(datnod(2,nl)) 
zl=abs(datnod(3,nl)) 
x2=abs (datnodCl ,n2) ) 
y2=abs (datnod(2 , n2) ) 
z2=abs (datnod( 3 , n2) ) 
xm=amaxl(xl,x2) 
ym=amaxl(yl,y2) 
zm=araaxl(zl ,z2) 
ix=0 
iy=0 
iz=0 

if (xm . le . edged) ix=igndp ( 1 ) 
if  (ym .  le .  edged)  iy=igndp(2) 
if (zm . le . edged) iz=igndp(3) 
npe  c=- (  aminO ( ix , 0 ) +amin0 ( iy , 0 ) +amin0 ( iz , 0 ) ) 
npmc=amaxO ( ix , 0 ) +amax0 ( iy , 0 ) +amax0 ( iz , 0 ) 
if (npec . ge . 1 ) then 

c  case  edge  is  attached  to  at  least  pec 
nunknb=nunknb+l 
nconn(3, ie)=nconn(3,ie)+l 
endif 

100  continue 
endif 

write(3 ,29) 

29  format (/14x, 'edge-vertex  connection  list'/) 
do  40  i=l,nedges 

write(3,331)i,nconn(l,i) ,nconn(2, i) ,nconn(3,i) 

40  continue 

331  format (3x, 'edge' ,i5, '  goes  from  vertex',i5,'  to  vertex',i5, 

$'  mult=',i3) 
return 
end 

subrout ine  prntbnd (nconn , nbound , ist art , i ,nedges , nf aces , nbodys ) 

c  this  subroutine  prints  the  edges  and  the  vertices  of  each  face. 

c  input : 

c  nconn  has  the  vertices  and  the  multiplicity  factor  for  each  edge. 

c  nbound  has  the  edges  for  each  face. 

c  ist  art  has  the  beginning  faces  for  each  body. 
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c  i  is  the  present  body, 
c  nedges  is  the  total  number  of  edges, 

c  nfaces  is  the  total  number  of  faces, 

c  nbodys  is  the  total  number  of  bodys. 

integer  nconn( 3, nedges) ,nbonnd(3, nfaces) ,istart (nbodys+1) 
integer  nverts(10000,3) 
open (20 , f ile=  '  node 1 . m  ' ) 
open ( 2 1 , f ile=  * node2 . m 1 ) 
open(22,f ile='node3.m') 
c  loop  through  the  faces  of  this  body. 

do  10  ilO=istart (i) ,istart (i+l)-l 

call  facvtx(nconn, nedges, nbound(l,il0) ,nbound(2, ilO) , 

>nbound(3,il0) ,nvl,nv2,nv3) 

write(3,98)il0,nbound(i,il0) ,nbound(2,il0) ,nbound(3,il0) ,nvl,nv2 
>,nv3 

c  write  face  number  and  vertices  to  "facedat" 
c  write(9,*)il0,nvi,nv2,nv3 
write (20, 100)  nvl 
write (21, 100)  nv2 
write (22, 100)  nv3 
nverts(il0, l)=nvl 
nverts(il0,2)=nv2 
nverts (ilO , 3) =nv3 
10  continue 

98  format (lx, 'face* ,i5, '  has  edges', 3i5,'  with  vertices ' ,3i5) 

100  format (i6) 
return 
end 

Subroutines  GEOM,  FACETCK,  and  AXB  are  used  here.  Listings  appear  in  Ap¬ 
pendix  A. 
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APPENDIX  G:  MATLAB  GEOMETRY  VIEWER  (PLTPATCH) 


PLTPATCH  uses  the  data  files  produced  by  BLDMAT  to  plot  the  “inpatch”  ge¬ 
ometry.  PLTPACH  has  flags  that  can  be  changed  by  the  user  to  choose  between  wire 
grid,  3D,  and  3D  with  hidden  lines. 

•/•  Program  to  read  "buildnS"  data  and  plot  it  in  MATLAB 
•/.  node  (x,y,z)  coordinates 
clear 
elf 

*/.  legend: 

V  icur=0:  plot  wire  grid 

*/.  icur=l:  surface  plot  (with  hidden  lines) 

Vi  icur=2:  surface  plot  with  current  intensities 

*/.  icur=3:  plot  currents  for  a  "deck  walk" 
icur=0; 
load  xpts.m 
load  ypts.m 
load  zpts.m 
nverts=length(xpts) ; 

Vi  node  connection  list 
load  endl.m 
load  end2.m 
nedges=G.ength(endl) ; 

Vi  face  connection  list 
load  nodel.m 
load  node 2. m 
load  node3.m 
nf aces=length(node3) ; 

Vi  load  vind  array  which  gives  the  three  nodes  for  each  triangle 
for  i=i:nfaces 

pts=[nodel(i)  node2(i)  node3(i)]  ; 
vind(i, : )=pts ; 
end 

Vi  Read  currents  if  intensities  are  to  be  plotted 
if  icur==2 

load  current s.ra 

currents=10*logl0(currents+l . e-10)  ; 
minc=min( currents) ; 
max  c  =iricLX  (  cur r  ent  s  )  ; 

current s  =  ( current  s -mine  )  /  ( maxc-minc )  ; 
end 

Vi  This  section  plots  a  mesh 
if  icur==0 
for  i=l:nfaces 

X= [xpts(vind(i, 1))  xpts(vind(i,2))  xpts(vind(i ,3) )  xpts(vind(i,l))] ; 

Y= [ypts(vind(i, 1))  ypts(vind(i,2))  3rpts(vind(i,3) )  ypts(vind(i,l))] ; 

Z= [zpts(vind(i, 1))  zpts(vind(i,2))  zpts(vind(i,3))  zpts(vind(i,l))] ; 
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plot3(X,Y,Z, ,m> ) 
if  i  ==  1 
axis  equal 
view(60,45) 
hold  on 
end 
end 
end 

'/.  This  section  plots  surfaces  (therefore  can  use  hidden) 
if  icur~=0 
'/.  length  64  map 
colormap(hot) 
cur  r  ent  s = curr  ent  s  *  64 ; 
it=0; 

for  i=l infaces 
it=it+i ; 

X= [xpts(vind(i,l))  xpts(vind(i, 1))  ;  xpts(vind(i,2))  xpts(vind(i,3))] ; 
Y=[ypts(vind(i,l))  ypts(vind(i,l))  ;  ypts(vind(i,2))  ypts(vind(i,3))] ; 
Z=[zpts(vind(i,l))  zpts(vind(i,l))  ;  zpts(vind(i,2))  zpts(vind(i,3))] ; 
*/•  Save  handles  for  each  triangle 
if  it==l 
axis  equal 
view(60,45) 
hold  on 
end 

if  icur==2 

C=[floor(currents(i))  f loor (currents (i))  ; 

f 1 o o r ( curr  ent  s ( i ) )  floor ( currents (i))] ; 
spit (it)=surf (X, Y,Z,C) ; 
end 

if  icur==l 

spit (it)=surf (X,Y,Z) ; 
set(splt(it)  ,  'facecolor' ,  'black')  ; 
end 

•/.  fill3(X,Y,Z,C); 

*/.  patch(x,y,z) 

set  (spit  (it)  ,  'edgecolor' , 'white')  ; 
end 
end 

7.  label  nodes  if  desired 
ilabn=l; 
if  ilabn==0 
for  i=l inverts 

text (xpts(i) ,ypts(i) ,zpts(i) ,num2str(i)) 
end 
end 

y,  label  edges 
ilabe=0; 
if  ilabe==0 
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for  i=l:nedges 

xav= ( xpt s ( endl ( i ) ) +xpt  s ( end2 ( i ) ) ) /2 . ; 
yav=(ypts (endl (i) )+ypts (end2(i) ) )/2 . ; 
zav=(zpts(endl(i))+zpts(end2(i)))/2. ; 
t  ext ( xav , yav , zav , num2  s tr ( i ) ) 
end 
end 
end 

axis  square 
xlabel( 'x') 
ylabel( *y * ) 
zlabelC'z1) 
hold  off 

delx=max(xpts)-min(xpts) ; 
dely=max(ypts)-min(ypts) ; 
delz=max(zpts)-rain(zpts) ; 
del=max(  [delx  dely  delz]  )  ; 

axis([min(xpts) ,min(xpts)+del,min(ypts),min(ypts)+del, . . . 
min(zpts) ,min(zpts)+del] ) 
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